Based on basic code for calculating and mapping climatic novelty using the sigma dissimilarity metric. * written by Colin Mahony * Supp Material of “A closer look at novel climates: new methods and insights at continental to landscape scales”

What this code does

Compare 1800 analog (A) to today surface (B), with today’s surface variability (C)

Compare today surface analog (A) to 2100-RCP8.5 (B), with today’s surface variability (C)

Compare today surface analog (A) to 2100-RCP4.5 (B), with today’s surface variability (C)

This code also compares the converse (e.g. the current climate as B to the future climate as A) to assess disappearing climates

This code uses the annnual climate normals and monthly ICV for the novelty calculation

##################################
#### System setup ####
##################################
## Specify location of data ####
setwd("/Users/lotterhos/Documents/GitHub/OceanClimateNovelty") 
#setwd("~/Desktop/PostDoc/OceanClimateNovelty/")

source("src/0-Novelty_Oceans_Functions.R")
## Warning: package 'raster' was built under R version 3.6.2
## Loading required package: sp
## Warning: package 'sp' was built under R version 3.6.2
## Loading required package: ade4
## Loading required package: adehabitatMA
## Registered S3 methods overwritten by 'adehabitatMA':
##   method                       from
##   print.SpatialPixelsDataFrame sp  
##   print.SpatialPixels          sp
## 
## Attaching package: 'adehabitatMA'
## The following object is masked from 'package:raster':
## 
##     buffer
## Loading required package: CircStats
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following objects are masked from 'package:raster':
## 
##     area, select
## Loading required package: boot
## 
## Attaching package: 'data.table'
## The following object is masked from 'package:raster':
## 
##     shift
## ── Attaching packages ────────────────────── tidyverse 1.2.1 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.3
## ✓ tidyr   1.0.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0
## Warning: package 'ggplot2' was built under R version 3.6.2
## ── Conflicts ───────────────────────── tidyverse_conflicts() ──
## x dplyr::between()   masks data.table::between()
## x tidyr::extract()   masks raster::extract()
## x dplyr::filter()    masks stats::filter()
## x dplyr::first()     masks data.table::first()
## x dplyr::id()        masks adehabitatLT::id()
## x dplyr::lag()       masks stats::lag()
## x dplyr::last()      masks data.table::last()
## x dplyr::select()    masks MASS::select(), raster::select()
## x purrr::transpose() masks data.table::transpose()
## Loading required package: spam
## Loading required package: dotCall64
## Loading required package: grid
## Spam version 2.4-0 (2019-11-01) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
## 
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
## 
##     backsolve, forwardsolve
## Loading required package: maps
## 
## Attaching package: 'maps'
## The following object is masked from 'package:purrr':
## 
##     map
## See https://github.com/NCAR/Fields for
##  an extensive vignette, other supplements and source code
## rgdal: version: 1.4-7, (SVN revision 845)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
##  Linking to sp version: 1.3-1
## Warning: replacing previous import 'sf::st_make_valid' by
## 'lwgeom::st_make_valid' when loading 'tmap'
## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo
## Checking rgeos availability: TRUE
## Warning: package 'sf' was built under R version 3.6.2
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
## 
## Attaching package: 'fasterize'
## The following object is masked from 'package:graphics':
## 
##     plot
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
## The following objects are masked from 'package:data.table':
## 
##     dcast, melt
##################################
#### Read in the input data ####
##################################

path <- "/Users/lotterhos/Google Drive/_LotterhosLab2021/LotLab_projects/0-grants_or_active/2018-Unfunded-OceanClimateNovelty/Data/Lotterhos/"

path2<- "/Users/lotterhos/Google Drive/_LotterhosLab2021/LotLab_projects/0-grants_or_active/2018-Unfunded-OceanClimateNovelty/Data/"

dat1 <- fread(paste0(path,"Katie_Temp_Arag_1800_2000.txt"), sep = ",")
dat4.5 <- fread(paste0(path,"Katie_Temp_Arag_2070_2100_RCP45.txt"), sep = ",")
dat8.5 <- fread(paste0(path,"Katie_Temp_Arag_2070_2100_RCP85.txt"), sep = ",")

ICV <- fread(paste0(path2, "ICV_Temp_Arag_pH_1960_2020.txt"))

results_path <- "/Users/lotterhos/Google Drive/_LotterhosLab2021/LotLab_projects/0-grants_or_active/2018-Unfunded-OceanClimateNovelty/Results/"
results_dir <- paste0(results_path, "HemisphereRestricted/")

head(dat1)
##    No  Lon   Lat Year Month      SST   Arag     pH
## 1:  1 20.5 -69.5 1800     1 -0.05589 1.8167 8.0967
## 2:  1 20.5 -69.5 1800     2 -0.49722 2.0572 8.1706
## 3:  1 20.5 -69.5 1800     3 -0.88001 1.9321 8.1444
## 4:  1 20.5 -69.5 1800     4 -1.24800 1.8865 8.1396
## 5:  1 20.5 -69.5 1800     5 -2.23610 1.8170 8.1375
## 6:  1 20.5 -69.5 1800     6 -0.21900 1.9615 8.1438
head(dat8.5)
##    No  Lon   Lat Year Month       SST    Arag     pH
## 1:  1 20.5 -69.5 2070     1  0.503480 0.95276 7.7819
## 2:  1 20.5 -69.5 2070     2  0.062159 1.00710 7.8161
## 3:  1 20.5 -69.5 2070     3 -0.320630 0.96601 7.8035
## 4:  1 20.5 -69.5 2070     4 -0.688590 0.94601 7.8008
## 5:  1 20.5 -69.5 2070     5 -1.676700 0.90633 7.7981
## 6:  1 20.5 -69.5 2070     6  0.340380 0.98976 7.8048
head(dat4.5)
##    No  Lon   Lat Year Month      SST   Arag     pH
## 1:  1 20.5 -69.5 2070     1  0.58164 1.1518 7.8685
## 2:  1 20.5 -69.5 2070     2  0.14032 1.2362 7.9107
## 3:  1 20.5 -69.5 2070     3 -0.24248 1.1804 7.8955
## 4:  1 20.5 -69.5 2070     4 -0.61043 1.1555 7.8924
## 5:  1 20.5 -69.5 2070     5 -1.59860 1.1084 7.8899
## 6:  1 20.5 -69.5 2070     6  0.41854 1.2072 7.8965
head(ICV)
##    No  Lon   Lat Year Month      SST   Arag     pH
## 1:  1 20.5 -69.5 1960     1 -0.20420 1.5702 8.0291
## 2:  1 20.5 -69.5 1960     2 -0.64552 1.7479 8.0921
## 3:  1 20.5 -69.5 1960     3 -1.02830 1.6505 8.0698
## 4:  1 20.5 -69.5 1960     4 -1.39630 1.6131 8.0657
## 5:  1 20.5 -69.5 1960     5 -2.38440 1.5515 8.0635
## 6:  1 20.5 -69.5 1960     6 -0.36730 1.6803 8.0699
unique(dat1$Year)
## [1] 1800 1810 1820 1830 1970 1980 1990 2000
unique(dat8.5$Year)
## [1] 2070 2080 2090 2100
unique(dat4.5$Year)
## [1] 2070 2080 2090 2100
unique(ICV$Year)
## [1] 1960 1970 1980 1990 2000 2010 2020
hist(c(dat1$Arag, dat8.5$Arag))

hist(log10(c(dat1$Arag, dat8.5$Arag)))

boxplot ICV

pdf(paste0(results_dir, "ICV_Lat.pdf"), width=9, height=7)
par(mfrow=c(3,1), mar=c(4,5,0.5,0.5), oma=c(2,0,0,0), las=0)

lat <- sort(unique(ICV$Lat))
length(lat)
## [1] 169
lat[1]
## [1] -78.5
lat[42]
## [1] -37.5
lat[80]
## [1] 0.5
lat[117]
## [1] 37.5
lat[158]
## [1] 78.5
ind <- c(1,42,80,117,158)
lat_names=rep(NA, 169)
lat_names[ind] <- lat[ind]
lat_names[80] <- 0 #close enough

boxplot(ICV$SST~ICV$Lat, col="green",  
        #border = rgb(0,0,0,0.3),  outpch=20, outcol= rgb(0,0,0,0.3), outbg=rgb(0,0,0,0.3), 
        ylab="", xlab="", names=lat_names, las=1, cex=1.2)
  mtext("SST", side=2, cex=1.5, line=3)

boxplot(ICV$Arag~ICV$Lat, col="green",
        ylab="", xlab="", names=lat_names, las=1, cex=1.2)
  mtext("Arag.", side=2, cex=1.5, line=3)

boxplot(ICV$pH ~ ICV$Lat, col="green",
        ylab="", xlab="", names=lat_names, las=1, cex=1.2)
  mtext("pH", side=2, cex=1.5, line=3)


mtext("Latitude", side=1, outer=TRUE, cex=2)
dev.off()
## quartz_off_screen 
##                 2
##################################
#### Log-transform Arag ####
##################################
dat1$Arag <- log10(dat1$Arag)
dat4.5$Arag <- log10(dat4.5$Arag)
dat8.5$Arag <- log10(dat8.5$Arag)
ICV$Arag <- log10(ICV$Arag)

#  Aragonite saturation is a ratio variable; it is limited at zero and 
# proportional changes are meaningful. The analysis uses raw values of 
# aragonite saturation state. Under raw scaling, the difference between 1 and 2 
# is equal to the difference between zero and one, which doesn’t reflect the 
# meaning of the variable. Instead, the difference between 1 and 2 should 
# be given the same significance as the difference between 0.5 and 1 (doubling vs halving). 
# Log-transforming aragonite saturation state will produce this more meaningful 
# proportional scaling.

# More generally on the same topic: The need for proportional scaling is the 
# reason why precipitation and other ratio variables (e.g. degree-days) are 
# typically log-transformed in climate space analysis. In theory, all variables 
# should be log-transformed to produce an environmental space where distances 
# are comparable for different locations; however, in practice temperature 
# doesn’t need to be log-transformed because it doesn’t vary across orders 
# of magnitude, and in our case pH is already a log-scaled variable.

#boxplot(dat8.5$pH ~ dat8.5$Year)
#boxplot(dat4.5$pH ~ dat4.5$Year)
# note that these two datasets are exactly the same for years < 2000
#--------------------------------
#### 40-year climate normals ####
#--------------------------------
# for example, 1930 represents from 1/1/1925 to 12/31/1934.  
dat_2000 <-   dat1 %>% filter(Year>1960 & Year<2010)
dim(dat_2000)
## [1] 1903776       8
# Below is a standardization code that can be deleted becuase the PCA
# does the standardization
# x_SST <- mean(dat_2000$SST)
# s_SST <- sd(dat_2000$SST)
# x_Arag <- mean(dat_2000$Arag)
# s_Arag <- sd(dat_2000$Arag)
# x_pH <- mean(dat_2000$pH)
# s_pH <- sd(dat_2000$pH)
# 
# head(dat_2000)  
#   dat_2000$SST <- (dat_2000$SST - x_SST)/s_SST
#   dat_2000$Arag <- (dat_2000$Arag - x_Arag)/s_Arag
#   dat_2000$pH <- (dat_2000$pH - x_pH)/s_pH
#   dat_2000 %>% summarise_each(sd)
# 
dat_1800 <- dat1 %>% filter(Year<1850)
dim(dat_1800)
## [1] 1903776       8
#   dat_1800$SST <- (dat_1800$SST - x_SST)/s_SST
#   dat_1800$Arag <- (dat_1800$Arag - x_Arag)/s_Arag
#   dat_1800$pH <- (dat_1800$pH - x_pH)/s_pH
#   
#   
dat_2100_8.5 <-   dat8.5 
#     dim(dat_2100_8.5) 
#     dat_2100_8.5$SST <- (dat_2100_8.5$SST - x_SST)/s_SST
#     dat_2100_8.5$Arag <- (dat_2100_8.5$Arag - x_Arag)/s_Arag
#     dat_2100_8.5$pH <- (dat_2100_8.5$pH - x_pH)/s_pH
#   
dat_2100_4.5 <-   dat4.5 
#   dim(dat_2100_4.5) 
#   dat_2100_4.5$SST <- (dat_2100_4.5$SST - x_SST)/s_SST
#   dat_2100_4.5$Arag <- (dat_2100_4.5$Arag - x_Arag)/s_Arag
#   dat_2100_4.5$pH <- (dat_2100_4.5$pH - x_pH)/s_pH

sum(!complete.cases(dat_1800))
## [1] 0
sum(!complete.cases(dat_2000))
## [1] 0
sum(!complete.cases(dat_2100_4.5))
## [1] 0
sum(!complete.cases(dat_2100_8.5))
## [1] 0
sum(!complete.cases(ICV))
## [1] 0
# Calculate climate normals
# Returns the mean climate for each station at each time
norm_1800 <- calculate_normals_annual(dat_1800)
norm_2000 <- calculate_normals_annual(dat_2000)
norm_2100_8.5 <- calculate_normals_annual(dat_2100_8.5)
norm_2100_4.5 <- calculate_normals_annual(dat_2100_4.5)

head(norm_1800)
##   No        SST      Arag       pH
## 1  1 -1.2720933 0.2661633 8.129290
## 2  2 -1.1802993 0.2581453 8.136665
## 3  3 -1.0816947 0.2502643 8.141752
## 4  4 -0.9757862 0.2425956 8.144575
## 5  5 -0.8618450 0.2352653 8.145356
## 6  6 -0.7388042 0.2284502 8.144419
dim(norm_1800)
## [1] 39662     4
dim(norm_2000)
## [1] 39662     4
dim(norm_2100_8.5)
## [1] 39662     4
dim(norm_2100_4.5)
## [1] 39662     4
#--------------------------------  
### data frame to link station number to Lat Long ####
#--------------------------------
head(dat8.5)
##    No  Lon   Lat Year Month       SST         Arag     pH
## 1:  1 20.5 -69.5 2070     1  0.503480 -0.021016484 7.7819
## 2:  1 20.5 -69.5 2070     2  0.062159  0.003072596 7.8161
## 3:  1 20.5 -69.5 2070     3 -0.320630 -0.015018378 7.8035
## 4:  1 20.5 -69.5 2070     4 -0.688590 -0.024104273 7.8008
## 5:  1 20.5 -69.5 2070     5 -1.676700 -0.042713644 7.7981
## 6:  1 20.5 -69.5 2070     6  0.340380 -0.004470102 7.8048
stations <- unique(dat8.5$No)
stationInfo = data.frame(stations=stations, lat=NA, long=NA)
unik <- which(!duplicated(dat8.5$No))
head(unik)
## [1]   1  49  97 145 193 241
stationInfo$lat <- dat8.5$Lat[unik]
stationInfo$long <- dat8.5$Lon[unik]
names(stationInfo)[1] <- "No"
head(stationInfo)
##   No   lat long
## 1  1 -69.5 20.5
## 2  2 -68.5 20.5
## 3  3 -67.5 20.5
## 4  4 -66.5 20.5
## 5  5 -65.5 20.5
## 6  6 -64.5 20.5
which(!complete.cases(stationInfo))
## integer(0)
#-------------------------------- 
### ICV matrix
head(ICV)  
##    No  Lon   Lat Year Month      SST      Arag     pH
## 1:  1 20.5 -69.5 1960     1 -0.20420 0.1959550 8.0291
## 2:  1 20.5 -69.5 1960     2 -0.64552 0.2425166 8.0921
## 3:  1 20.5 -69.5 1960     3 -1.02830 0.2176155 8.0698
## 4:  1 20.5 -69.5 1960     4 -1.39630 0.2076613 8.0657
## 5:  1 20.5 -69.5 1960     5 -2.38440 0.1907518 8.0635
## 6:  1 20.5 -69.5 1960     6 -0.36730 0.2253868 8.0699
C <- ICV %>% select(No, SST, Arag, pH)
head(C)
##    No      SST      Arag     pH
## 1:  1 -0.20420 0.1959550 8.0291
## 2:  1 -0.64552 0.2425166 8.0921
## 3:  1 -1.02830 0.2176155 8.0698
## 4:  1 -1.39630 0.2076613 8.0657
## 5:  1 -2.38440 0.1907518 8.0635
## 6:  1 -0.36730 0.2253868 8.0699
#--------------------------------   
#---------------
### Hemisphere Indexes ###
# N. hemisphere
N_hem <- which(stationInfo$lat>=0)
S_hem <- which(stationInfo$lat<0)
#---------------

1800 analog to today (Novelty)

#--------------------------------  
### 1800 analog to today (Novelty) ####
# What are today's novel climates compared to 1800?
# Novel climates are identified by comparing the 
# future climate realization (B) for each gridpoint 
# to the past (A) climate realizations (Williams)
# RESTRICT BASELINE TO HEMISPHERE OF OBSERVATION
#--------------------------------
# A is baseline, the base that a station from B is compared to
# If A is past and B is future, then degree of novelty
# If A is future and B is past, then degree of disappearance
# B will be subset to a particular station for analysis
# C is the data frame used to calculate the ICV, will also be subset to a particular station for analysis


A <- norm_1800
# 1800 climate normals
# Serves as the Baseline
head(A)
##   No        SST      Arag       pH
## 1  1 -1.2720933 0.2661633 8.129290
## 2  2 -1.1802993 0.2581453 8.136665
## 3  3 -1.0816947 0.2502643 8.141752
## 4  4 -0.9757862 0.2425956 8.144575
## 5  5 -0.8618450 0.2352653 8.145356
## 6  6 -0.7388042 0.2284502 8.144419
dim(A)
## [1] 39662     4
# subset the 2000 data to the same stations as the 1800 data
B <- norm_2000
# 2000 climate normals 
head(B)
##   No        SST      Arag       pH
## 1  1 -1.4479254 0.1799026 8.035921
## 2  2 -1.3696480 0.1741309 8.045892
## 3  3 -1.2845583 0.1683092 8.053425
## 4  4 -1.1916660 0.1624748 8.058488
## 5  5 -1.0895625 0.1567233 8.061231
## 6  6 -0.9764333 0.1512065 8.062000
dim(B)
## [1] 39662     4
# same columns as A

# sanity check to make sure stations in right order
identical(A$No, B$No) # should be true
## [1] TRUE
# check list 18906,18952, 29736, 29789, 8193, 29319,
# calc_sigma_D(A, B, C,  18906, "_1800_today")
# calc_sigma_D(A, B, C,  18952, "_1800_today")
# calc_sigma_D(A, B, C,  29736, "_1800_today")
# calc_sigma_D(A, B, C,  29789, "_1800_today")
# calc_sigma_D(A, B, C,  8193, "_1800_today")
# calc_sigma_D(A, B, C,  29319, "_1800_today")


#NN.sigma_1800_today <- loop_sigma_D(A, B, C, "_1800_today") # uses whole ocean as a baseline
NN.sigma_1800_today_N <- loop_sigma_D(A[N_hem,], B[N_hem,], C, "_1800_today")
## [1]   100   317 17413
## [1]   200   610 17413
## [1]   300   934 17413
## [1]   400  1247 17413
## [1]   500  1683 17413
## [1]   600  2053 17413
## [1]   700  2356 17413
## [1]   800  2660 17413
## [1]   900  2971 17413
## [1]  1000  3213 17413
## [1]  1100  3452 17413
## [1]  1200  3759 17413
## [1]  1300  3996 17413
## [1]  1400  4236 17413
## [1]  1500  4540 17413
## [1]  1600  4847 17413
## [1]  1700  5220 17413
## [1]  1800  5530 17413
## [1]  1900  5830 17413
## [1]  2000  6129 17413
## [1]  2100  6427 17413
## [1]  2200  6726 17413
## [1]  2300  7099 17413
## [1]  2400  7694 17413
## [1]  2500  8267 17413
## [1]  2600  8629 17413
## [1]  2700  8887 17413
## [1]  2800  9069 17413
## [1]  2900  9286 17413
## [1]  3000  9454 17413
## [1]  3100  9622 17413
## [1]  3200  9790 17413
## [1]  3300  9954 17413
## [1]  3400 10119 17413
## [1]  3500 10287 17413
## [1]  3600 10458 17413
## [1]  3700 10592 17413
## [1]  3800 10764 17413
## [1]  3900 10945 17413
## [1]  4000 11095 17413
## [1]  4100 11309 17413
## [1]  4200 11469 17413
## [1]  4300 11702 17413
## [1]  4400 11871 17413
## [1]  4500 12109 17413
## [1]  4600 12278 17413
## [1]  4700 12450 17413
## [1]  4800 12692 17413
## [1]  4900 12863 17413
## [1]  5000 13035 17413
## [1]  5100 13280 17413
## [1]  5200 13453 17413
## [1]  5300 13627 17413
## [1]  5400 13878 17413
## [1]  5500 14056 17413
## [1]  5600 14235 17413
## [1]  5700 14414 17413
## [1]  5800 14670 17413
## [1]  5900 14848 17413
## [1]  6000 15026 17413
## [1]  6100 15204 17413
## [1]  6200 15384 17413
## [1]  6300 15641 17413
## [1]  6400 15819 17413
## [1]  6500 15998 17413
## [1]  6600 16176 17413
## [1]  6700 16356 17413
## [1]  6800 16534 17413
## [1]  6900 16791 17413
## [1]  7000 16970 17413
## [1]  7100 17149 17413
## [1]  7200 17328 17413
## [1]  7300 17506 17413
## [1]  7400 17684 17413
## [1]  7500 17862 17413
## [1]  7600 18119 17413
## [1]  7700 18297 17413
## [1]  7800 18476 17413
## [1]  7900 18655 17413
## [1]  8000 18912 17413
## [1]  8100 19090 17413
## [1]  8200 19269 17413
## [1]  8300 19526 17413
## [1]  8400 19704 17413
## [1]  8500 19881 17413
## [1]  8600 20058 17413
## [1]  8700 20312 17413
## [1]  8800 20488 17413
## [1]  8900 20664 17413
## [1]  9000 20916 17413
## [1]  9100 21091 17413
## [1]  9200 21266 17413
## [1]  9300 21441 17413
## [1]  9400 21691 17413
## [1]  9500 21866 17413
## [1]  9600 22041 17413
## [1]  9700 22218 17413
## [1]  9800 22468 17413
## [1]  9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413
## [1] 12300 29551 17413
## [1] 12400 29721 17413
## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413
NN.sigma_1800_today_S <- loop_sigma_D(A[S_hem,], B[S_hem,], C, "_1800_today")
## [1]   100   140 22249
## [1]   200   298 22249
## [1]   300   459 22249
## [1]   400   622 22249
## [1]   500   765 22249
## [1]   600   909 22249
## [1]   700  1052 22249
## [1]   800  1195 22249
## [1]   900  1345 22249
## [1]  1000  1468 22249
## [1]  1100  1619 22249
## [1]  1200  1763 22249
## [1]  1300  1923 22249
## [1]  1400  2093 22249
## [1]  1500  2228 22249
## [1]  1600  2401 22249
## [1]  1700  2539 22249
## [1]  1800  2718 22249
## [1]  1900  2858 22249
## [1]  2000  3040 22249
## [1]  2100  3182 22249
## [1]  2200  3372 22249
## [1]  2300  3518 22249
## [1]  2400  3710 22249
## [1]  2500  3856 22249
## [1]  2600  4047 22249
## [1]  2700  4190 22249
## [1]  2800  4368 22249
## [1]  2900  4506 22249
## [1]  3000  4678 22249
## [1]  3100  4810 22249
## [1]  3200  4967 22249
## [1]  3300  5092 22249
## [1]  3400  5240 22249
## [1]  3500  5367 22249
## [1]  3600  5531 22249
## [1]  3700  5666 22249
## [1]  3800  5834 22249
## [1]  3900  5969 22249
## [1]  4000  6141 22249
## [1]  4100  6278 22249
## [1]  4200  6451 22249
## [1]  4300  6586 22249
## [1]  4400  6753 22249
## [1]  4500  6882 22249
## [1]  4600  7031 22249
## [1]  4700  7149 22249
## [1]  4800  7275 22249
## [1]  4900  7388 22249
## [1]  5000  7514 22249
## [1]  5100  7640 22249
## [1]  5200  7753 22249
## [1]  5300  7882 22249
## [1]  5400  8010 22249
## [1]  5500  8139 22249
## [1]  5600  8256 22249
## [1]  5700  8404 22249
## [1]  5800  8538 22249
## [1]  5900  8680 22249
## [1]  6000  8864 22249
## [1]  6100  9073 22249
## [1]  6200  9259 22249
## [1]  6300  9501 22249
## [1]  6400  9758 22249
## [1]  6500 10019 22249
## [1]  6600 10336 22249
## [1]  6700 10614 22249
## [1]  6800 10899 22249
## [1]  6900 11128 22249
## [1]  7000 11363 22249
## [1]  7100 11532 22249
## [1]  7200 11769 22249
## [1]  7300 11938 22249
## [1]  7400 12178 22249
## [1]  7500 12350 22249
## [1]  7600 12523 22249
## [1]  7700 12772 22249
## [1]  7800 12948 22249
## [1]  7900 13203 22249
## [1]  8000 13382 22249
## [1]  8100 13561 22249
## [1]  8200 13824 22249
## [1]  8300 14012 22249
## [1]  8400 14279 22249
## [1]  8500 14462 22249
## [1]  8600 14649 22249
## [1]  8700 14915 22249
## [1]  8800 15102 22249
## [1]  8900 15289 22249
## [1]  9000 15556 22249
## [1]  9100 15740 22249
## [1]  9200 15926 22249
## [1]  9300 16112 22249
## [1]  9400 16386 22249
## [1]  9500 16573 22249
## [1]  9600 16760 22249
## [1]  9700 16947 22249
## [1]  9800 17224 22249
## [1]  9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249
NN.sigma_1800_today <- rbind(NN.sigma_1800_today_N, NN.sigma_1800_today_S)
head(NN.sigma_1800_today)
##   No NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1 36         0.033725924                 30751          0.23355016
## 2 37         0.015960753                 30751          0.16010043
## 3 38         0.010016616                 31111          0.12668130
## 4 39         0.001600598                 14563          0.05055509
## 5 40         0.018324556                  1958          0.17162744
## 6 41         0.004808260                 17177          0.08767895
##   numPCs_1800_today
## 1                 2
## 2                 2
## 3                 2
## 4                 2
## 5                 2
## 6                 2
hist(NN.sigma_1800_today$numPCs_1800_today)

final_dat <- full_join(stationInfo, NN.sigma_1800_today,  by="No")
head(final_dat)
##   No   lat long NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1  1 -69.5 20.5         0.005340623                 40066           0.0924152
## 2  2 -68.5 20.5         0.007501424                 40009           0.1095737
## 3  3 -67.5 20.5         0.009629160                 40067           0.1241975
## 4  4 -66.5 20.5         0.015490226                 40009           0.1577081
## 5  5 -65.5 20.5         0.024201364                 40066           0.1974682
## 6  6 -64.5 20.5         0.049860481                 40066           0.2848815
##   numPCs_1800_today
## 1                 2
## 2                 2
## 3                 2
## 4                 2
## 5                 2
## 6                 2
final_dat[which(is.infinite(final_dat$NN.sigma_1800_today)),]
## [1] No                    lat                   long                 
## [4] NN.sigma_1800_today   NN.station_1800_today NN.Mdist_1800_today  
## [7] numPCs_1800_today    
## <0 rows> (or 0-length row.names)
final_dat[which(!complete.cases(final_dat)),]
## [1] No                    lat                   long                 
## [4] NN.sigma_1800_today   NN.station_1800_today NN.Mdist_1800_today  
## [7] numPCs_1800_today    
## <0 rows> (or 0-length row.names)
#NN.sigma_1800_today[which(is.infinite(NN.sigma_1800_today))] <- NA

identical(nrow(final_dat), nrow(stationInfo)) 
## [1] TRUE
# should be true, same number as stations as stationInfo

# Visualize
world <- map_data("world2")
Plot_nonInt(final_dat$lat, final_dat$long, 
            final_dat$NN.sigma_1800_today, world, "sigD novel 2000")

# The NN.sigma_1800_today represents the novelty of today's
# climate compared to 1800
# NN.dist_1800_today represents the Mahalanobis distance from the climate
# in today to it's nearest neighbor in all 1800 data
# No--> NN.station_1800_today represents the station that was queried 2000 (B) 
# and it's nearest data in the global reference NN.station_1800_today (1800, A)
# SEE CODE BELOW FOR ADDITION OF NN SITES

Today analog to 1800 Disappearance

#--------------------------------  
### Today analog to 1800 Disappearance####
# What are 1800 disappearing climates?
# disappearing climates are identified by comparing
# each past gridpoint (A) to all future climate realizations (B)
#--------------------------------
calc_sigma_D(A, B, C, 1)
##      NN.sigma NN.station   NN.Mdist num_PCs
## 1 0.005073206      22827 0.09006697       2
calc_sigma_D(B, A, C, 1)
##       NN.sigma NN.station   NN.Mdist num_PCs
## 1 0.0005525579      38202 0.02969763       2
B <- norm_1800
A <- norm_2000
identical(A$No, B$No) # should be true
## [1] TRUE
identical(sort(A$No), A$No) # should be true
## [1] TRUE
# check list 18906,18952, 29736, 29789, 8193, 29319,
# calc_sigma_D(A, B, C,  18906, "_today_1800")
# calc_sigma_D(A, B, C,  18952, "_today_1800")
# calc_sigma_D(A, B, C,  29736, "_today_1800")
# calc_sigma_D(A, B, C,  29789, "_today_1800")
# calc_sigma_D(A, B, C,  8193, "_today_1800")
# calc_sigma_D(A, B, C,  29319, "_today_1800")


NN.sigma_today_1800_N <- loop_sigma_D(A[N_hem,], B[N_hem,], C, "_today_1800")
## [1]   100   317 17413
## [1]   200   610 17413
## [1]   300   934 17413
## [1]   400  1247 17413
## [1]   500  1683 17413
## [1]   600  2053 17413
## [1]   700  2356 17413
## [1]   800  2660 17413
## [1]   900  2971 17413
## [1]  1000  3213 17413
## [1]  1100  3452 17413
## [1]  1200  3759 17413
## [1]  1300  3996 17413
## [1]  1400  4236 17413
## [1]  1500  4540 17413
## [1]  1600  4847 17413
## [1]  1700  5220 17413
## [1]  1800  5530 17413
## [1]  1900  5830 17413
## [1]  2000  6129 17413
## [1]  2100  6427 17413
## [1]  2200  6726 17413
## [1]  2300  7099 17413
## [1]  2400  7694 17413
## [1]  2500  8267 17413
## [1]  2600  8629 17413
## [1]  2700  8887 17413
## [1]  2800  9069 17413
## [1]  2900  9286 17413
## [1]  3000  9454 17413
## [1]  3100  9622 17413
## [1]  3200  9790 17413
## [1]  3300  9954 17413
## [1]  3400 10119 17413
## [1]  3500 10287 17413
## [1]  3600 10458 17413
## [1]  3700 10592 17413
## [1]  3800 10764 17413
## [1]  3900 10945 17413
## [1]  4000 11095 17413
## [1]  4100 11309 17413
## [1]  4200 11469 17413
## [1]  4300 11702 17413
## [1]  4400 11871 17413
## [1]  4500 12109 17413
## [1]  4600 12278 17413
## [1]  4700 12450 17413
## [1]  4800 12692 17413
## [1]  4900 12863 17413
## [1]  5000 13035 17413
## [1]  5100 13280 17413
## [1]  5200 13453 17413
## [1]  5300 13627 17413
## [1]  5400 13878 17413
## [1]  5500 14056 17413
## [1]  5600 14235 17413
## [1]  5700 14414 17413
## [1]  5800 14670 17413
## [1]  5900 14848 17413
## [1]  6000 15026 17413
## [1]  6100 15204 17413
## [1]  6200 15384 17413
## [1]  6300 15641 17413
## [1]  6400 15819 17413
## [1]  6500 15998 17413
## [1]  6600 16176 17413
## [1]  6700 16356 17413
## [1]  6800 16534 17413
## [1]  6900 16791 17413
## [1]  7000 16970 17413
## [1]  7100 17149 17413
## [1]  7200 17328 17413
## [1]  7300 17506 17413
## [1]  7400 17684 17413
## [1]  7500 17862 17413
## [1]  7600 18119 17413
## [1]  7700 18297 17413
## [1]  7800 18476 17413
## [1]  7900 18655 17413
## [1]  8000 18912 17413
## [1]  8100 19090 17413
## [1]  8200 19269 17413
## [1]  8300 19526 17413
## [1]  8400 19704 17413
## [1]  8500 19881 17413
## [1]  8600 20058 17413
## [1]  8700 20312 17413
## [1]  8800 20488 17413
## [1]  8900 20664 17413
## [1]  9000 20916 17413
## [1]  9100 21091 17413
## [1]  9200 21266 17413
## [1]  9300 21441 17413
## [1]  9400 21691 17413
## [1]  9500 21866 17413
## [1]  9600 22041 17413
## [1]  9700 22218 17413
## [1]  9800 22468 17413
## [1]  9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413
## [1] 12300 29551 17413
## [1] 12400 29721 17413
## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413
NN.sigma_today_1800_S <- loop_sigma_D(A[S_hem,], B[S_hem,], C, "_today_1800")
## [1]   100   140 22249
## [1]   200   298 22249
## [1]   300   459 22249
## [1]   400   622 22249
## [1]   500   765 22249
## [1]   600   909 22249
## [1]   700  1052 22249
## [1]   800  1195 22249
## [1]   900  1345 22249
## [1]  1000  1468 22249
## [1]  1100  1619 22249
## [1]  1200  1763 22249
## [1]  1300  1923 22249
## [1]  1400  2093 22249
## [1]  1500  2228 22249
## [1]  1600  2401 22249
## [1]  1700  2539 22249
## [1]  1800  2718 22249
## [1]  1900  2858 22249
## [1]  2000  3040 22249
## [1]  2100  3182 22249
## [1]  2200  3372 22249
## [1]  2300  3518 22249
## [1]  2400  3710 22249
## [1]  2500  3856 22249
## [1]  2600  4047 22249
## [1]  2700  4190 22249
## [1]  2800  4368 22249
## [1]  2900  4506 22249
## [1]  3000  4678 22249
## [1]  3100  4810 22249
## [1]  3200  4967 22249
## [1]  3300  5092 22249
## [1]  3400  5240 22249
## [1]  3500  5367 22249
## [1]  3600  5531 22249
## [1]  3700  5666 22249
## [1]  3800  5834 22249
## [1]  3900  5969 22249
## [1]  4000  6141 22249
## [1]  4100  6278 22249
## [1]  4200  6451 22249
## [1]  4300  6586 22249
## [1]  4400  6753 22249
## [1]  4500  6882 22249
## [1]  4600  7031 22249
## [1]  4700  7149 22249
## [1]  4800  7275 22249
## [1]  4900  7388 22249
## [1]  5000  7514 22249
## [1]  5100  7640 22249
## [1]  5200  7753 22249
## [1]  5300  7882 22249
## [1]  5400  8010 22249
## [1]  5500  8139 22249
## [1]  5600  8256 22249
## [1]  5700  8404 22249
## [1]  5800  8538 22249
## [1]  5900  8680 22249
## [1]  6000  8864 22249
## [1]  6100  9073 22249
## [1]  6200  9259 22249
## [1]  6300  9501 22249
## [1]  6400  9758 22249
## [1]  6500 10019 22249
## [1]  6600 10336 22249
## [1]  6700 10614 22249
## [1]  6800 10899 22249
## [1]  6900 11128 22249
## [1]  7000 11363 22249
## [1]  7100 11532 22249
## [1]  7200 11769 22249
## [1]  7300 11938 22249
## [1]  7400 12178 22249
## [1]  7500 12350 22249
## [1]  7600 12523 22249
## [1]  7700 12772 22249
## [1]  7800 12948 22249
## [1]  7900 13203 22249
## [1]  8000 13382 22249
## [1]  8100 13561 22249
## [1]  8200 13824 22249
## [1]  8300 14012 22249
## [1]  8400 14279 22249
## [1]  8500 14462 22249
## [1]  8600 14649 22249
## [1]  8700 14915 22249
## [1]  8800 15102 22249
## [1]  8900 15289 22249
## [1]  9000 15556 22249
## [1]  9100 15740 22249
## [1]  9200 15926 22249
## [1]  9300 16112 22249
## [1]  9400 16386 22249
## [1]  9500 16573 22249
## [1]  9600 16760 22249
## [1]  9700 16947 22249
## [1]  9800 17224 22249
## [1]  9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249
NN.sigma_today_1800 <- rbind(NN.sigma_today_1800_N,
                             NN.sigma_today_1800_S)

head(NN.sigma_today_1800)
##   No NN.sigma_today_1800 NN.station_today_1800 NN.Mdist_today_1800
## 1 36         0.836709022                 31931           1.3486467
## 2 37         0.894133247                   794           1.4077485
## 3 38         0.661532267                   861           1.1633920
## 4 39         0.265489647                 39602           0.6854521
## 5 40         0.006257156                 39603           0.1000496
## 6 41         0.036266868                 39426           0.2423104
##   numPCs_today_1800
## 1                 2
## 2                 2
## 3                 2
## 4                 2
## 5                 2
## 6                 2
final_dat2 <- full_join(final_dat, NN.sigma_today_1800)
## Joining, by = "No"
head(final_dat2)
##   No   lat long NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1  1 -69.5 20.5         0.005340623                 40066           0.0924152
## 2  2 -68.5 20.5         0.007501424                 40009           0.1095737
## 3  3 -67.5 20.5         0.009629160                 40067           0.1241975
## 4  4 -66.5 20.5         0.015490226                 40009           0.1577081
## 5  5 -65.5 20.5         0.024201364                 40066           0.1974682
## 6  6 -64.5 20.5         0.049860481                 40066           0.2848815
##   numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## 1                 2         0.019575390                 26505
## 2                 2         0.020981196                 29981
## 3                 2         0.014089519                 26082
## 4                 2         0.006465863                 26082
## 5                 2         0.011645698                 25749
## 6                 2         0.005112661                 25973
##   NN.Mdist_today_1800 numPCs_today_1800
## 1          0.17743251                 2
## 2          0.18374455                 2
## 3          0.15036689                 2
## 4          0.10170870                 2
## 5          0.13663930                 2
## 6          0.09041724                 2
dim(final_dat2)
## [1] 39662    11
Plot_nonInt(final_dat2$lat, final_dat2$long, 
            final_dat2$NN.sigma_today_1800, world, "sigD dis 2000")

ggplot(final_dat2) + geom_point(aes(x=NN.sigma_1800_today, 
                                    y=NN.sigma_today_1800,
                                    color=lat), alpha=0.5) +
  scale_color_gradient2(low="red", high="blue", mid="grey") +
  theme_classic() + xlab("Novel climates today from 1800") +
  ylab("Disappearing 1800 climates")

# The NN.sigma_1800_disappear represents the novelty of 1800
# compared to today (e.g disappearing climates)
# NN.dist_1800_disappear represents the Mahalanobis distance from the climate
# in 1800 to it's nearest neighbor in all of today's data
# No--> NN.station_1800_disappear represents the the station that was queried 1800 (B) 
# and it's nearest data in the global reference NN.station_1800_disappear (2000, A)

Add the info for lat/long of nearest neighbors (1800 to today)

stationInfo2 <- stationInfo
names(stationInfo2) <- paste0(names(stationInfo2),"_1800_today")
head(stationInfo2)
##   No_1800_today lat_1800_today long_1800_today
## 1             1          -69.5            20.5
## 2             2          -68.5            20.5
## 3             3          -67.5            20.5
## 4             4          -66.5            20.5
## 5             5          -65.5            20.5
## 6             6          -64.5            20.5
dim(stationInfo)
## [1] 39662     3
names(stationInfo2)[1] <-"NN.station_1800_today"
final_dat3 <- left_join(final_dat2, stationInfo2)
## Joining, by = "NN.station_1800_today"
dim(final_dat3)
## [1] 39662    13
head(final_dat3)
##   No   lat long NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1  1 -69.5 20.5         0.005340623                 40066           0.0924152
## 2  2 -68.5 20.5         0.007501424                 40009           0.1095737
## 3  3 -67.5 20.5         0.009629160                 40067           0.1241975
## 4  4 -66.5 20.5         0.015490226                 40009           0.1577081
## 5  5 -65.5 20.5         0.024201364                 40066           0.1974682
## 6  6 -64.5 20.5         0.049860481                 40066           0.2848815
##   numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## 1                 2         0.019575390                 26505
## 2                 2         0.020981196                 29981
## 3                 2         0.014089519                 26082
## 4                 2         0.006465863                 26082
## 5                 2         0.011645698                 25749
## 6                 2         0.005112661                 25973
##   NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today long_1800_today
## 1          0.17743251                 2          -69.5           379.5
## 2          0.18374455                 2          -69.5           378.5
## 3          0.15036689                 2          -68.5           379.5
## 4          0.10170870                 2          -69.5           378.5
## 5          0.13663930                 2          -69.5           379.5
## 6          0.09041724                 2          -69.5           379.5
tail(final_dat3)
##          No  lat  long NN.sigma_1800_today NN.station_1800_today
## 39657 40115 84.5 379.5        0.0033782225                 37019
## 39658 40116 85.5 379.5        0.0041649559                 37769
## 39659 40117 86.5 379.5        0.0001037297                  3229
## 39660 40118 87.5 379.5        0.0047175381                  7551
## 39661 40119 88.5 379.5        0.0377153031                  7551
## 39662 40120 89.5 379.5        0.0646634655                  7551
##       NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 39657          0.07347198                 2            1.728258
## 39658          0.08159263                 2            1.931682
## 39659          0.01286606                 2            2.170640
## 39660          0.08684629                 2            2.350300
## 39661          0.24717272                 2            2.424026
## 39662          0.32537576                 2            2.455741
##       NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800
## 39657                   613            2.226041                 2
## 39658                   613            2.420731                 2
## 39659                   613            2.648753                 2
## 39660                   554            2.819972                 2
## 39661                  3223            2.890216                 2
## 39662                  3336            2.920432                 2
##       lat_1800_today long_1800_today
## 39657           86.5           347.5
## 39658           86.5           353.5
## 39659           88.5            59.5
## 39660           88.5           103.5
## 39661           88.5           103.5
## 39662           88.5           103.5
# Should both be 0
sum(!complete.cases(final_dat2))
## [1] 0
sum(!complete.cases(final_dat3))
## [1] 0
head(final_dat3[!complete.cases(final_dat3),])
##  [1] No                    lat                   long                 
##  [4] NN.sigma_1800_today   NN.station_1800_today NN.Mdist_1800_today  
##  [7] numPCs_1800_today     NN.sigma_today_1800   NN.station_today_1800
## [10] NN.Mdist_today_1800   numPCs_today_1800     lat_1800_today       
## [13] long_1800_today      
## <0 rows> (or 0-length row.names)
summary(final_dat3, NN.station_today_1800)
##        No             lat               long       NN.sigma_1800_today
##  Min.   :    1   Min.   :-78.500   Min.   : 20.5   Min.   :0.0000002  
##  1st Qu.:10100   1st Qu.:-43.500   1st Qu.:138.5   1st Qu.:0.0037303  
##  Median :20110   Median : -9.500   Median :207.5   Median :0.0778035  
##  Mean   :20095   Mean   : -1.731   Mean   :208.8   Mean   :0.4287141  
##  3rd Qu.:30087   3rd Qu.: 34.500   3rd Qu.:295.5   3rd Qu.:0.6417823  
##  Max.   :40120   Max.   : 89.500   Max.   :379.5   Max.   :3.1374782  
##  NN.station_1800_today NN.Mdist_1800_today numPCs_1800_today
##  Min.   :   44         Min.   :0.000572    Min.   :2.000    
##  1st Qu.:20665         1st Qu.:0.077212    1st Qu.:2.000    
##  Median :25736         Median :0.357832    Median :2.000    
##  Mean   :24568         Mean   :0.705606    Mean   :2.002    
##  3rd Qu.:29342         3rd Qu.:1.142396    3rd Qu.:2.000    
##  Max.   :40074         Max.   :3.611830    Max.   :3.000    
##  NN.sigma_today_1800 NN.station_today_1800 NN.Mdist_today_1800
##  Min.   :0.00000     Min.   :   39         Min.   :0.000314   
##  1st Qu.:0.03089     1st Qu.:17767         1st Qu.:0.224749   
##  Median :0.83348     Median :30424         Median :1.345494   
##  Mean   :0.92699     Mean   :24144         Mean   :1.296300   
##  3rd Qu.:1.60850     3rd Qu.:31698         3rd Qu.:2.112351   
##  Max.   :3.44851     Max.   :40119         Max.   :3.868075   
##  numPCs_today_1800 lat_1800_today    long_1800_today
##  Min.   :2.000     Min.   :-77.500   Min.   : 20.5  
##  1st Qu.:2.000     1st Qu.:-43.500   1st Qu.:210.5  
##  Median :2.000     Median : -2.500   Median :246.5  
##  Mean   :2.002     Mean   : -2.377   Mean   :244.0  
##  3rd Qu.:2.000     3rd Qu.: 25.500   3rd Qu.:284.5  
##  Max.   :3.000     Max.   : 89.500   Max.   :379.5
final_dat3$latshift_1800_today <- NA
condN <- which(final_dat3$lat>0)
final_dat3$latshift_1800_today[condN] <- final_dat3$lat[condN] - 
  final_dat3$lat_1800_today[condN]
condS <- which(final_dat3$lat<0)
final_dat3$latshift[condS] <- -1*(final_dat3$lat[condS] - 
                              final_dat3$lat_1800_today[condS])
hist(final_dat3$latshift_1800_today)

#stationInfo3 <- stationInfo
#names(stationInfo3) <- paste0(names(stationInfo3),"_today_1800")
#names(stationInfo3)[1] <- "NN.station_today_1800"
#final_dat4 <- left_join(final_dat3, stationInfo3)
#dim(final_dat4)
#head(final_dat4)

# cond <- which(final_dat4$NN.sigma_1800_today>5 |  final_dat4$NN.sigma_1800_disappear>5)
# length(cond)
# final_dat4[cond,]
ggplot(final_dat3) + geom_point(aes(y=lat, 
                                    x=lat_1800_today,
                                    color=NN.sigma_1800_today), alpha=0.2) +
  scale_color_gradient2(low="red", high="blue", mid="grey") +
  theme_classic() + ylab("Latitude in 2000") +
  xlab("Latitude of nearest neighbor in 1800") + geom_abline(intercept=0,slope=1) 

# This plot makes more sense in the way we typically think about it.
# If we consider a location in 2000, where was the climate it came from in 1800?


write.csv(final_dat3, paste0(results_dir, "1800base_2000_novelty_NShemrestrict.csv"))
final_dat4 <- final_dat3
# ggplot(final_dat4) + geom_point(aes(x=lat_1800_disappear,
#                                     y=lat,
#                                     color=NN.sigma_1800_disappear), alpha=0.2) +
#   scale_color_gradient2(low="red", high="blue", mid="grey") +
#   theme_classic() + ylab("Latitude in 1800") +
#   xlab("Latitude of nearest neighbor in 2000") + geom_abline(intercept=0,slope=1) +
#   geom_abline(intercept=0,slope=-1)
# this is weird

#n <- final_dat4 %>% filter(lat_1800_disappear > 0 & lat_1800_disappear < 10)
#Plot_nonInt(n$lat_1800_disappear, n$long_1800_disappear, 
#            n$NN.sigma_1800_disappear, world, "sigma dis.")
#hist(final_dat4$lat_1800_disappear, breaks=seq(-90,90, by=1))
#hist(final_dat4$lat_1800_today, breaks=seq(-90,90, by=1))
#hist(final_dat4$lat, breaks=seq(-90,90, by=1))

Today analog to 2100 RCP 8.5 (Novelty)

identical(norm_2000$No, norm_2100_8.5$No)
## [1] TRUE
A1 <- norm_2000
# 1970-2000 climate normals
head(A1)
##   No        SST      Arag       pH
## 1  1 -1.4479254 0.1799026 8.035921
## 2  2 -1.3696480 0.1741309 8.045892
## 3  3 -1.2845583 0.1683092 8.053425
## 4  4 -1.1916660 0.1624748 8.058488
## 5  5 -1.0895625 0.1567233 8.061231
## 6  6 -0.9764333 0.1512065 8.062000
dim(A1)
## [1] 39662     4
B1 <- norm_2100_8.5
# 2070-2100
head(B1)
##   No         SST        Arag       pH
## 1  1 -0.55466154 -0.09210795 7.726612
## 2  2 -0.45603873 -0.09383454 7.740267
## 3  3 -0.34939250 -0.09577437 7.751387
## 4  4 -0.23345842 -0.09793483 7.759896
## 5  5 -0.10654381 -0.10027507 7.765894
## 6  6  0.03343595 -0.10269243 7.769621
dim(B1)
## [1] 39662     4
# same columns as A

# sanity check to make sure stations in right order
identical(A1$No, B1$No) # should be true
## [1] TRUE
identical(A1$No, sort(A1$No)) # should also be true
## [1] TRUE
# check list 18906,18952, 29736, 29789, 8193, 29319,
# calc_sigma_D(A1, B1, C,  18906, "_today_2100_8.5")
# calc_sigma_D(A1, B1, C,  18952, "_today_2100_8.5")
# calc_sigma_D(A1, B1, C,  29736, "_today_2100_8.5")
# calc_sigma_D(A1, B1, C,  29789, "_today_2100_8.5")
# calc_sigma_D(A1, B1, C,  8193, "_today_2100_8.5")
# calc_sigma_D(A1, B1, C,  29319, "_today_2100_8.5")

# C is the same defined above
NN.sigma_today_2100_8.5_N <- loop_sigma_D(A1[N_hem,], B1[N_hem,], C, "_today_2100_8.5")
## [1]   100   317 17413
## [1]   200   610 17413
## [1]   300   934 17413
## [1]   400  1247 17413
## [1]   500  1683 17413
## [1]   600  2053 17413
## [1]   700  2356 17413
## [1]   800  2660 17413
## [1]   900  2971 17413
## [1]  1000  3213 17413
## [1]  1100  3452 17413
## [1]  1200  3759 17413
## [1]  1300  3996 17413
## [1]  1400  4236 17413
## [1]  1500  4540 17413
## [1]  1600  4847 17413
## [1]  1700  5220 17413
## [1]  1800  5530 17413
## [1]  1900  5830 17413
## [1]  2000  6129 17413
## [1]  2100  6427 17413
## [1]  2200  6726 17413
## [1]  2300  7099 17413
## [1]  2400  7694 17413
## [1]  2500  8267 17413
## [1]  2600  8629 17413
## [1]  2700  8887 17413
## [1]  2800  9069 17413
## [1]  2900  9286 17413
## [1]  3000  9454 17413
## [1]  3100  9622 17413
## [1]  3200  9790 17413
## [1]  3300  9954 17413
## [1]  3400 10119 17413
## [1]  3500 10287 17413
## [1]  3600 10458 17413
## [1]  3700 10592 17413
## [1]  3800 10764 17413
## [1]  3900 10945 17413
## [1]  4000 11095 17413
## [1]  4100 11309 17413
## [1]  4200 11469 17413
## [1]  4300 11702 17413
## [1]  4400 11871 17413
## [1]  4500 12109 17413
## [1]  4600 12278 17413
## [1]  4700 12450 17413
## [1]  4800 12692 17413
## [1]  4900 12863 17413
## [1]  5000 13035 17413
## [1]  5100 13280 17413
## [1]  5200 13453 17413
## [1]  5300 13627 17413
## [1]  5400 13878 17413
## [1]  5500 14056 17413
## [1]  5600 14235 17413
## [1]  5700 14414 17413
## [1]  5800 14670 17413
## [1]  5900 14848 17413
## [1]  6000 15026 17413
## [1]  6100 15204 17413
## [1]  6200 15384 17413
## [1]  6300 15641 17413
## [1]  6400 15819 17413
## [1]  6500 15998 17413
## [1]  6600 16176 17413
## [1]  6700 16356 17413
## [1]  6800 16534 17413
## [1]  6900 16791 17413
## [1]  7000 16970 17413
## [1]  7100 17149 17413
## [1]  7200 17328 17413
## [1]  7300 17506 17413
## [1]  7400 17684 17413
## [1]  7500 17862 17413
## [1]  7600 18119 17413
## [1]  7700 18297 17413
## [1]  7800 18476 17413
## [1]  7900 18655 17413
## [1]  8000 18912 17413
## [1]  8100 19090 17413
## [1]  8200 19269 17413
## [1]  8300 19526 17413
## [1]  8400 19704 17413
## [1]  8500 19881 17413
## [1]  8600 20058 17413
## [1]  8700 20312 17413
## [1]  8800 20488 17413
## [1]  8900 20664 17413
## [1]  9000 20916 17413
## [1]  9100 21091 17413
## [1]  9200 21266 17413
## [1]  9300 21441 17413
## [1]  9400 21691 17413
## [1]  9500 21866 17413
## [1]  9600 22041 17413
## [1]  9700 22218 17413
## [1]  9800 22468 17413
## [1]  9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413
## [1] 12300 29551 17413
## [1] 12400 29721 17413
## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413
# C is the same defined above
NN.sigma_today_2100_8.5_S <- loop_sigma_D(A1[S_hem,], B1[S_hem,], C, "_today_2100_8.5")
## [1]   100   140 22249
## [1]   200   298 22249
## [1]   300   459 22249
## [1]   400   622 22249
## [1]   500   765 22249
## [1]   600   909 22249
## [1]   700  1052 22249
## [1]   800  1195 22249
## [1]   900  1345 22249
## [1]  1000  1468 22249
## [1]  1100  1619 22249
## [1]  1200  1763 22249
## [1]  1300  1923 22249
## [1]  1400  2093 22249
## [1]  1500  2228 22249
## [1]  1600  2401 22249
## [1]  1700  2539 22249
## [1]  1800  2718 22249
## [1]  1900  2858 22249
## [1]  2000  3040 22249
## [1]  2100  3182 22249
## [1]  2200  3372 22249
## [1]  2300  3518 22249
## [1]  2400  3710 22249
## [1]  2500  3856 22249
## [1]  2600  4047 22249
## [1]  2700  4190 22249
## [1]  2800  4368 22249
## [1]  2900  4506 22249
## [1]  3000  4678 22249
## [1]  3100  4810 22249
## [1]  3200  4967 22249
## [1]  3300  5092 22249
## [1]  3400  5240 22249
## [1]  3500  5367 22249
## [1]  3600  5531 22249
## [1]  3700  5666 22249
## [1]  3800  5834 22249
## [1]  3900  5969 22249
## [1]  4000  6141 22249
## [1]  4100  6278 22249
## [1]  4200  6451 22249
## [1]  4300  6586 22249
## [1]  4400  6753 22249
## [1]  4500  6882 22249
## [1]  4600  7031 22249
## [1]  4700  7149 22249
## [1]  4800  7275 22249
## [1]  4900  7388 22249
## [1]  5000  7514 22249
## [1]  5100  7640 22249
## [1]  5200  7753 22249
## [1]  5300  7882 22249
## [1]  5400  8010 22249
## [1]  5500  8139 22249
## [1]  5600  8256 22249
## [1]  5700  8404 22249
## [1]  5800  8538 22249
## [1]  5900  8680 22249
## [1]  6000  8864 22249
## [1]  6100  9073 22249
## [1]  6200  9259 22249
## [1]  6300  9501 22249
## [1]  6400  9758 22249
## [1]  6500 10019 22249
## [1]  6600 10336 22249
## [1]  6700 10614 22249
## [1]  6800 10899 22249
## [1]  6900 11128 22249
## [1]  7000 11363 22249
## [1]  7100 11532 22249
## [1]  7200 11769 22249
## [1]  7300 11938 22249
## [1]  7400 12178 22249
## [1]  7500 12350 22249
## [1]  7600 12523 22249
## [1]  7700 12772 22249
## [1]  7800 12948 22249
## [1]  7900 13203 22249
## [1]  8000 13382 22249
## [1]  8100 13561 22249
## [1]  8200 13824 22249
## [1]  8300 14012 22249
## [1]  8400 14279 22249
## [1]  8500 14462 22249
## [1]  8600 14649 22249
## [1]  8700 14915 22249
## [1]  8800 15102 22249
## [1]  8900 15289 22249
## [1]  9000 15556 22249
## [1]  9100 15740 22249
## [1]  9200 15926 22249
## [1]  9300 16112 22249
## [1]  9400 16386 22249
## [1]  9500 16573 22249
## [1]  9600 16760 22249
## [1]  9700 16947 22249
## [1]  9800 17224 22249
## [1]  9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249
NN.sigma_today_2100_8.5 <- rbind(NN.sigma_today_2100_8.5_N, 
                                 NN.sigma_today_2100_8.5_S)


which(is.infinite(NN.sigma_today_2100_8.5$NN.sigma_today_2100_8.5))
## integer(0)
which(is.na(NN.sigma_today_2100_8.5$NN.sigma_today_2100_8.5))
## integer(0)
final_dat5 <- full_join(NN.sigma_today_2100_8.5, final_dat4, by="No")
head(final_dat5)
##   No NN.sigma_today_2100_8.5 NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5
## 1 36            0.3501681236                     29638              0.79989070
## 2 37            0.1043489359                     29695              0.41656957
## 3 38            0.0015820951                     29638              0.05026184
## 4 39            0.0006816715                     28978              0.03298615
## 5 40            0.1188184652                     28440              0.44577407
## 6 41            0.5463419465                     28254              1.03579199
##   numPCs_today_2100_8.5  lat long NN.sigma_1800_today NN.station_1800_today
## 1                     2 70.5 20.5         0.033725924                 30751
## 2                     2 71.5 20.5         0.015960753                 30751
## 3                     2 72.5 20.5         0.010016616                 31111
## 4                     2 73.5 20.5         0.001600598                 14563
## 5                     2 74.5 20.5         0.018324556                  1958
## 6                     2 75.5 20.5         0.004808260                 17177
##   NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 1          0.23355016                 2         0.836709022
## 2          0.16010043                 2         0.894133247
## 3          0.12668130                 2         0.661532267
## 4          0.05055509                 2         0.265489647
## 5          0.17162744                 2         0.006257156
## 6          0.08767895                 2         0.036266868
##   NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## 1                 31931           1.3486467                 2           44.5
## 2                   794           1.4077485                 2           44.5
## 3                   861           1.1633920                 2           82.5
## 4                 39602           0.6854521                 2           53.5
## 5                 39603           0.1000496                 2           69.5
## 6                 39426           0.2423104                 2           61.5
##   long_1800_today latshift_1800_today latshift
## 1           302.5                  26       26
## 2           302.5                  27       27
## 3           305.5                 -10      -10
## 4           172.5                  20       20
## 5            47.5                   5        5
## 6           188.5                  14       14
dim(final_dat5)
## [1] 39662    19
stationInfo5 <- stationInfo
names(stationInfo5) <- paste0(names(stationInfo),"_today_2100_8.5")
names(stationInfo5)[1] <- "NN.station_today_2100_8.5"
head(stationInfo5)
##   NN.station_today_2100_8.5 lat_today_2100_8.5 long_today_2100_8.5
## 1                         1              -69.5                20.5
## 2                         2              -68.5                20.5
## 3                         3              -67.5                20.5
## 4                         4              -66.5                20.5
## 5                         5              -65.5                20.5
## 6                         6              -64.5                20.5
final_dat6 <- left_join(final_dat5, stationInfo5)
## Joining, by = "NN.station_today_2100_8.5"
dim(final_dat6)
## [1] 39662    21
head(final_dat6)
##   No NN.sigma_today_2100_8.5 NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5
## 1 36            0.3501681236                     29638              0.79989070
## 2 37            0.1043489359                     29695              0.41656957
## 3 38            0.0015820951                     29638              0.05026184
## 4 39            0.0006816715                     28978              0.03298615
## 5 40            0.1188184652                     28440              0.44577407
## 6 41            0.5463419465                     28254              1.03579199
##   numPCs_today_2100_8.5  lat long NN.sigma_1800_today NN.station_1800_today
## 1                     2 70.5 20.5         0.033725924                 30751
## 2                     2 71.5 20.5         0.015960753                 30751
## 3                     2 72.5 20.5         0.010016616                 31111
## 4                     2 73.5 20.5         0.001600598                 14563
## 5                     2 74.5 20.5         0.018324556                  1958
## 6                     2 75.5 20.5         0.004808260                 17177
##   NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 1          0.23355016                 2         0.836709022
## 2          0.16010043                 2         0.894133247
## 3          0.12668130                 2         0.661532267
## 4          0.05055509                 2         0.265489647
## 5          0.17162744                 2         0.006257156
## 6          0.08767895                 2         0.036266868
##   NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## 1                 31931           1.3486467                 2           44.5
## 2                   794           1.4077485                 2           44.5
## 3                   861           1.1633920                 2           82.5
## 4                 39602           0.6854521                 2           53.5
## 5                 39603           0.1000496                 2           69.5
## 6                 39426           0.2423104                 2           61.5
##   long_1800_today latshift_1800_today latshift lat_today_2100_8.5
## 1           302.5                  26       26               89.5
## 2           302.5                  27       27               89.5
## 3           305.5                 -10      -10               89.5
## 4           172.5                  20       20               85.5
## 5            47.5                   5        5               82.5
## 6           188.5                  14       14               82.5
##   long_today_2100_8.5
## 1               287.5
## 2               288.5
## 3               287.5
## 4               279.5
## 5               273.5
## 6               271.5
# Visualize
Plot_nonInt(final_dat6$lat, final_dat6$long, 
            final_dat6$NN.sigma_today_2100_8.5, world, "sigma nov 2100 8.5")

2100 RCP 8.5 to today, what are today’s climates that will disappear in 2100 RCP 8.5?

#--------------------------------  
### 2100 RCP 8.5 to today, what are today's climates that ####
### will disappear in 2100 RCP 8.5?
#--------------------------------
# check list 18906,18952, 29736, 29789, 8193, 29319,
# calc_sigma_D( B1, A1, C,  18906, "_2100_8.5_today")
# calc_sigma_D( B1, A1, C,  18952, "_2100_8.5_today")
# calc_sigma_D( B1, A1, C,  29736, "_2100_8.5_today")
# calc_sigma_D( B1, A1, C,  29789, "_2100_8.5_today")
# calc_sigma_D( B1, A1, C,  8193, "_2100_8.5_today")
# calc_sigma_D( B1, A1, C,  29319, "_2100_8.5_today")

NN.sigma_2100_8.5_today_N <- loop_sigma_D( B1[N_hem,], A1[N_hem,], C, "_2100_8.5_today")
## [1]   100   317 17413
## [1]   200   610 17413
## [1]   300   934 17413
## [1]   400  1247 17413
## [1]   500  1683 17413
## [1]   600  2053 17413
## [1]   700  2356 17413
## [1]   800  2660 17413
## [1]   900  2971 17413
## [1]  1000  3213 17413
## [1]  1100  3452 17413
## [1]  1200  3759 17413
## [1]  1300  3996 17413
## [1]  1400  4236 17413
## [1]  1500  4540 17413
## [1]  1600  4847 17413
## [1]  1700  5220 17413
## [1]  1800  5530 17413
## [1]  1900  5830 17413
## [1]  2000  6129 17413
## [1]  2100  6427 17413
## [1]  2200  6726 17413
## [1]  2300  7099 17413
## [1]  2400  7694 17413
## [1]  2500  8267 17413
## [1]  2600  8629 17413
## [1]  2700  8887 17413
## [1]  2800  9069 17413
## [1]  2900  9286 17413
## [1]  3000  9454 17413
## [1]  3100  9622 17413
## [1]  3200  9790 17413
## [1]  3300  9954 17413
## [1]  3400 10119 17413
## [1]  3500 10287 17413
## [1]  3600 10458 17413
## [1]  3700 10592 17413
## [1]  3800 10764 17413
## [1]  3900 10945 17413
## [1]  4000 11095 17413
## [1]  4100 11309 17413
## [1]  4200 11469 17413
## [1]  4300 11702 17413
## [1]  4400 11871 17413
## [1]  4500 12109 17413
## [1]  4600 12278 17413
## [1]  4700 12450 17413
## [1]  4800 12692 17413
## [1]  4900 12863 17413
## [1]  5000 13035 17413
## [1]  5100 13280 17413
## [1]  5200 13453 17413
## [1]  5300 13627 17413
## [1]  5400 13878 17413
## [1]  5500 14056 17413
## [1]  5600 14235 17413
## [1]  5700 14414 17413
## [1]  5800 14670 17413
## [1]  5900 14848 17413
## [1]  6000 15026 17413
## [1]  6100 15204 17413
## [1]  6200 15384 17413
## [1]  6300 15641 17413
## [1]  6400 15819 17413
## [1]  6500 15998 17413
## [1]  6600 16176 17413
## [1]  6700 16356 17413
## [1]  6800 16534 17413
## [1]  6900 16791 17413
## [1]  7000 16970 17413
## [1]  7100 17149 17413
## [1]  7200 17328 17413
## [1]  7300 17506 17413
## [1]  7400 17684 17413
## [1]  7500 17862 17413
## [1]  7600 18119 17413
## [1]  7700 18297 17413
## [1]  7800 18476 17413
## [1]  7900 18655 17413
## [1]  8000 18912 17413
## [1]  8100 19090 17413
## [1]  8200 19269 17413
## [1]  8300 19526 17413
## [1]  8400 19704 17413
## [1]  8500 19881 17413
## [1]  8600 20058 17413
## [1]  8700 20312 17413
## [1]  8800 20488 17413
## [1]  8900 20664 17413
## [1]  9000 20916 17413
## [1]  9100 21091 17413
## [1]  9200 21266 17413
## [1]  9300 21441 17413
## [1]  9400 21691 17413
## [1]  9500 21866 17413
## [1]  9600 22041 17413
## [1]  9700 22218 17413
## [1]  9800 22468 17413
## [1]  9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413
## [1] 12300 29551 17413
## [1] 12400 29721 17413
## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413
NN.sigma_2100_8.5_today_S <- loop_sigma_D( B1[S_hem,], A1[S_hem,], C, "_2100_8.5_today")
## [1]   100   140 22249
## [1]   200   298 22249
## [1]   300   459 22249
## [1]   400   622 22249
## [1]   500   765 22249
## [1]   600   909 22249
## [1]   700  1052 22249
## [1]   800  1195 22249
## [1]   900  1345 22249
## [1]  1000  1468 22249
## [1]  1100  1619 22249
## [1]  1200  1763 22249
## [1]  1300  1923 22249
## [1]  1400  2093 22249
## [1]  1500  2228 22249
## [1]  1600  2401 22249
## [1]  1700  2539 22249
## [1]  1800  2718 22249
## [1]  1900  2858 22249
## [1]  2000  3040 22249
## [1]  2100  3182 22249
## [1]  2200  3372 22249
## [1]  2300  3518 22249
## [1]  2400  3710 22249
## [1]  2500  3856 22249
## [1]  2600  4047 22249
## [1]  2700  4190 22249
## [1]  2800  4368 22249
## [1]  2900  4506 22249
## [1]  3000  4678 22249
## [1]  3100  4810 22249
## [1]  3200  4967 22249
## [1]  3300  5092 22249
## [1]  3400  5240 22249
## [1]  3500  5367 22249
## [1]  3600  5531 22249
## [1]  3700  5666 22249
## [1]  3800  5834 22249
## [1]  3900  5969 22249
## [1]  4000  6141 22249
## [1]  4100  6278 22249
## [1]  4200  6451 22249
## [1]  4300  6586 22249
## [1]  4400  6753 22249
## [1]  4500  6882 22249
## [1]  4600  7031 22249
## [1]  4700  7149 22249
## [1]  4800  7275 22249
## [1]  4900  7388 22249
## [1]  5000  7514 22249
## [1]  5100  7640 22249
## [1]  5200  7753 22249
## [1]  5300  7882 22249
## [1]  5400  8010 22249
## [1]  5500  8139 22249
## [1]  5600  8256 22249
## [1]  5700  8404 22249
## [1]  5800  8538 22249
## [1]  5900  8680 22249
## [1]  6000  8864 22249
## [1]  6100  9073 22249
## [1]  6200  9259 22249
## [1]  6300  9501 22249
## [1]  6400  9758 22249
## [1]  6500 10019 22249
## [1]  6600 10336 22249
## [1]  6700 10614 22249
## [1]  6800 10899 22249
## [1]  6900 11128 22249
## [1]  7000 11363 22249
## [1]  7100 11532 22249
## [1]  7200 11769 22249
## [1]  7300 11938 22249
## [1]  7400 12178 22249
## [1]  7500 12350 22249
## [1]  7600 12523 22249
## [1]  7700 12772 22249
## [1]  7800 12948 22249
## [1]  7900 13203 22249
## [1]  8000 13382 22249
## [1]  8100 13561 22249
## [1]  8200 13824 22249
## [1]  8300 14012 22249
## [1]  8400 14279 22249
## [1]  8500 14462 22249
## [1]  8600 14649 22249
## [1]  8700 14915 22249
## [1]  8800 15102 22249
## [1]  8900 15289 22249
## [1]  9000 15556 22249
## [1]  9100 15740 22249
## [1]  9200 15926 22249
## [1]  9300 16112 22249
## [1]  9400 16386 22249
## [1]  9500 16573 22249
## [1]  9600 16760 22249
## [1]  9700 16947 22249
## [1]  9800 17224 22249
## [1]  9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249
NN.sigma_2100_8.5_today <- rbind(NN.sigma_2100_8.5_today_N, NN.sigma_2100_8.5_today_S)


final_dat7 <- full_join(NN.sigma_2100_8.5_today, final_dat6, by="No")
head(final_dat7)
##   No NN.sigma_2100_8.5_today NN.station_2100_8.5_today NN.Mdist_2100_8.5_today
## 1 36                8.292361                     40111               10.055106
## 2 37                8.292361                     40111               10.340310
## 3 38                8.292361                     40111               10.400029
## 4 39                8.292361                     40111               10.280321
## 5 40                8.292361                     40111               10.091702
## 6 41                8.292361                     40111                9.955604
##   numPCs_2100_8.5_today NN.sigma_today_2100_8.5 NN.station_today_2100_8.5
## 1                     2            0.3501681236                     29638
## 2                     2            0.1043489359                     29695
## 3                     2            0.0015820951                     29638
## 4                     2            0.0006816715                     28978
## 5                     2            0.1188184652                     28440
## 6                     2            0.5463419465                     28254
##   NN.Mdist_today_2100_8.5 numPCs_today_2100_8.5  lat long NN.sigma_1800_today
## 1              0.79989070                     2 70.5 20.5         0.033725924
## 2              0.41656957                     2 71.5 20.5         0.015960753
## 3              0.05026184                     2 72.5 20.5         0.010016616
## 4              0.03298615                     2 73.5 20.5         0.001600598
## 5              0.44577407                     2 74.5 20.5         0.018324556
## 6              1.03579199                     2 75.5 20.5         0.004808260
##   NN.station_1800_today NN.Mdist_1800_today numPCs_1800_today
## 1                 30751          0.23355016                 2
## 2                 30751          0.16010043                 2
## 3                 31111          0.12668130                 2
## 4                 14563          0.05055509                 2
## 5                  1958          0.17162744                 2
## 6                 17177          0.08767895                 2
##   NN.sigma_today_1800 NN.station_today_1800 NN.Mdist_today_1800
## 1         0.836709022                 31931           1.3486467
## 2         0.894133247                   794           1.4077485
## 3         0.661532267                   861           1.1633920
## 4         0.265489647                 39602           0.6854521
## 5         0.006257156                 39603           0.1000496
## 6         0.036266868                 39426           0.2423104
##   numPCs_today_1800 lat_1800_today long_1800_today latshift_1800_today latshift
## 1                 2           44.5           302.5                  26       26
## 2                 2           44.5           302.5                  27       27
## 3                 2           82.5           305.5                 -10      -10
## 4                 2           53.5           172.5                  20       20
## 5                 2           69.5            47.5                   5        5
## 6                 2           61.5           188.5                  14       14
##   lat_today_2100_8.5 long_today_2100_8.5
## 1               89.5               287.5
## 2               89.5               288.5
## 3               89.5               287.5
## 4               85.5               279.5
## 5               82.5               273.5
## 6               82.5               271.5
Plot_nonInt(final_dat7$lat, final_dat7$long, 
            final_dat7$NN.sigma_2100_8.5_today, world, "sigma dis. 2000 in 2100 8.5")

# 
# stationInfo7 <- stationInfo
# names(stationInfo7) <- paste0(names(stationInfo),"_2100_8.5_today")
# names(stationInfo7)[1] <- "NN.station_2100_8.5_today"
# head(stationInfo7)
# final_dat8 <- left_join(final_dat7, stationInfo7)
final_dat8 <- final_dat7

2100 RCP 4.5 against today (novelty)

#--------------------------------  
### 2100 RCP 4.5 against today (novelty) ####
#--------------------------------
identical(norm_2000$No, norm_2100_4.5$No)
## [1] TRUE
# A1 the same as above

B2 <- norm_2100_4.5
# 2070-2100
head(B2)
##   No        SST       Arag       pH
## 1  1 -0.8282522 0.04708031 7.879460
## 2  2 -0.7305017 0.04380220 7.891665
## 3  3 -0.6256332 0.04038376 7.901365
## 4  4 -0.5128039 0.03683358 7.908494
## 5  5 -0.3908187 0.03319882 7.913204
## 6  6 -0.2581154 0.02960423 7.915740
dim(B2)
## [1] 39662     4
# same columns as A

# sanity check to make sure stations in right order
identical(A1$No, B2$No) # should be true
## [1] TRUE
identical(A1$No, sort(A1$No)) # should also be true
## [1] TRUE
# C is the same defined above
NN.sigma_today_2100_4.5_N <- loop_sigma_D(A1[N_hem,], B2[N_hem,], C, "_today_2100_4.5")
## [1]   100   317 17413
## [1]   200   610 17413
## [1]   300   934 17413
## [1]   400  1247 17413
## [1]   500  1683 17413
## [1]   600  2053 17413
## [1]   700  2356 17413
## [1]   800  2660 17413
## [1]   900  2971 17413
## [1]  1000  3213 17413
## [1]  1100  3452 17413
## [1]  1200  3759 17413
## [1]  1300  3996 17413
## [1]  1400  4236 17413
## [1]  1500  4540 17413
## [1]  1600  4847 17413
## [1]  1700  5220 17413
## [1]  1800  5530 17413
## [1]  1900  5830 17413
## [1]  2000  6129 17413
## [1]  2100  6427 17413
## [1]  2200  6726 17413
## [1]  2300  7099 17413
## [1]  2400  7694 17413
## [1]  2500  8267 17413
## [1]  2600  8629 17413
## [1]  2700  8887 17413
## [1]  2800  9069 17413
## [1]  2900  9286 17413
## [1]  3000  9454 17413
## [1]  3100  9622 17413
## [1]  3200  9790 17413
## [1]  3300  9954 17413
## [1]  3400 10119 17413
## [1]  3500 10287 17413
## [1]  3600 10458 17413
## [1]  3700 10592 17413
## [1]  3800 10764 17413
## [1]  3900 10945 17413
## [1]  4000 11095 17413
## [1]  4100 11309 17413
## [1]  4200 11469 17413
## [1]  4300 11702 17413
## [1]  4400 11871 17413
## [1]  4500 12109 17413
## [1]  4600 12278 17413
## [1]  4700 12450 17413
## [1]  4800 12692 17413
## [1]  4900 12863 17413
## [1]  5000 13035 17413
## [1]  5100 13280 17413
## [1]  5200 13453 17413
## [1]  5300 13627 17413
## [1]  5400 13878 17413
## [1]  5500 14056 17413
## [1]  5600 14235 17413
## [1]  5700 14414 17413
## [1]  5800 14670 17413
## [1]  5900 14848 17413
## [1]  6000 15026 17413
## [1]  6100 15204 17413
## [1]  6200 15384 17413
## [1]  6300 15641 17413
## [1]  6400 15819 17413
## [1]  6500 15998 17413
## [1]  6600 16176 17413
## [1]  6700 16356 17413
## [1]  6800 16534 17413
## [1]  6900 16791 17413
## [1]  7000 16970 17413
## [1]  7100 17149 17413
## [1]  7200 17328 17413
## [1]  7300 17506 17413
## [1]  7400 17684 17413
## [1]  7500 17862 17413
## [1]  7600 18119 17413
## [1]  7700 18297 17413
## [1]  7800 18476 17413
## [1]  7900 18655 17413
## [1]  8000 18912 17413
## [1]  8100 19090 17413
## [1]  8200 19269 17413
## [1]  8300 19526 17413
## [1]  8400 19704 17413
## [1]  8500 19881 17413
## [1]  8600 20058 17413
## [1]  8700 20312 17413
## [1]  8800 20488 17413
## [1]  8900 20664 17413
## [1]  9000 20916 17413
## [1]  9100 21091 17413
## [1]  9200 21266 17413
## [1]  9300 21441 17413
## [1]  9400 21691 17413
## [1]  9500 21866 17413
## [1]  9600 22041 17413
## [1]  9700 22218 17413
## [1]  9800 22468 17413
## [1]  9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413
## [1] 12300 29551 17413
## [1] 12400 29721 17413
## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413
NN.sigma_today_2100_4.5_S <- loop_sigma_D(A1[S_hem,], B2[S_hem,], C, "_today_2100_4.5")
## [1]   100   140 22249
## [1]   200   298 22249
## [1]   300   459 22249
## [1]   400   622 22249
## [1]   500   765 22249
## [1]   600   909 22249
## [1]   700  1052 22249
## [1]   800  1195 22249
## [1]   900  1345 22249
## [1]  1000  1468 22249
## [1]  1100  1619 22249
## [1]  1200  1763 22249
## [1]  1300  1923 22249
## [1]  1400  2093 22249
## [1]  1500  2228 22249
## [1]  1600  2401 22249
## [1]  1700  2539 22249
## [1]  1800  2718 22249
## [1]  1900  2858 22249
## [1]  2000  3040 22249
## [1]  2100  3182 22249
## [1]  2200  3372 22249
## [1]  2300  3518 22249
## [1]  2400  3710 22249
## [1]  2500  3856 22249
## [1]  2600  4047 22249
## [1]  2700  4190 22249
## [1]  2800  4368 22249
## [1]  2900  4506 22249
## [1]  3000  4678 22249
## [1]  3100  4810 22249
## [1]  3200  4967 22249
## [1]  3300  5092 22249
## [1]  3400  5240 22249
## [1]  3500  5367 22249
## [1]  3600  5531 22249
## [1]  3700  5666 22249
## [1]  3800  5834 22249
## [1]  3900  5969 22249
## [1]  4000  6141 22249
## [1]  4100  6278 22249
## [1]  4200  6451 22249
## [1]  4300  6586 22249
## [1]  4400  6753 22249
## [1]  4500  6882 22249
## [1]  4600  7031 22249
## [1]  4700  7149 22249
## [1]  4800  7275 22249
## [1]  4900  7388 22249
## [1]  5000  7514 22249
## [1]  5100  7640 22249
## [1]  5200  7753 22249
## [1]  5300  7882 22249
## [1]  5400  8010 22249
## [1]  5500  8139 22249
## [1]  5600  8256 22249
## [1]  5700  8404 22249
## [1]  5800  8538 22249
## [1]  5900  8680 22249
## [1]  6000  8864 22249
## [1]  6100  9073 22249
## [1]  6200  9259 22249
## [1]  6300  9501 22249
## [1]  6400  9758 22249
## [1]  6500 10019 22249
## [1]  6600 10336 22249
## [1]  6700 10614 22249
## [1]  6800 10899 22249
## [1]  6900 11128 22249
## [1]  7000 11363 22249
## [1]  7100 11532 22249
## [1]  7200 11769 22249
## [1]  7300 11938 22249
## [1]  7400 12178 22249
## [1]  7500 12350 22249
## [1]  7600 12523 22249
## [1]  7700 12772 22249
## [1]  7800 12948 22249
## [1]  7900 13203 22249
## [1]  8000 13382 22249
## [1]  8100 13561 22249
## [1]  8200 13824 22249
## [1]  8300 14012 22249
## [1]  8400 14279 22249
## [1]  8500 14462 22249
## [1]  8600 14649 22249
## [1]  8700 14915 22249
## [1]  8800 15102 22249
## [1]  8900 15289 22249
## [1]  9000 15556 22249
## [1]  9100 15740 22249
## [1]  9200 15926 22249
## [1]  9300 16112 22249
## [1]  9400 16386 22249
## [1]  9500 16573 22249
## [1]  9600 16760 22249
## [1]  9700 16947 22249
## [1]  9800 17224 22249
## [1]  9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249
NN.sigma_today_2100_4.5 <- rbind(NN.sigma_today_2100_4.5_N, NN.sigma_today_2100_4.5_S)

which(is.infinite(NN.sigma_today_2100_4.5$NN.sigma_today_2100_4.5))
## integer(0)
which(is.na(NN.sigma_today_2100_4.5$NN.sigma_today_2100_4.5))
## integer(0)
final_dat9 <- full_join(NN.sigma_today_2100_4.5, final_dat8, by="No")
head(final_dat9)
##   No NN.sigma_today_2100_4.5 NN.station_today_2100_4.5 NN.Mdist_today_2100_4.5
## 1 36             0.260091777                     31342              0.67775158
## 2 37             0.150575943                     31338              0.50493740
## 3 38             0.005534165                     31456              0.09407847
## 4 39             0.139495959                     30347              0.48495881
## 5 40             0.011775140                      3437              0.13740012
## 6 41             0.098007656                      4227              0.40321232
##   numPCs_today_2100_4.5 NN.sigma_2100_8.5_today NN.station_2100_8.5_today
## 1                     2                8.292361                     40111
## 2                     2                8.292361                     40111
## 3                     2                8.292361                     40111
## 4                     2                8.292361                     40111
## 5                     2                8.292361                     40111
## 6                     2                8.292361                     40111
##   NN.Mdist_2100_8.5_today numPCs_2100_8.5_today NN.sigma_today_2100_8.5
## 1               10.055106                     2            0.3501681236
## 2               10.340310                     2            0.1043489359
## 3               10.400029                     2            0.0015820951
## 4               10.280321                     2            0.0006816715
## 5               10.091702                     2            0.1188184652
## 6                9.955604                     2            0.5463419465
##   NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5 numPCs_today_2100_8.5  lat
## 1                     29638              0.79989070                     2 70.5
## 2                     29695              0.41656957                     2 71.5
## 3                     29638              0.05026184                     2 72.5
## 4                     28978              0.03298615                     2 73.5
## 5                     28440              0.44577407                     2 74.5
## 6                     28254              1.03579199                     2 75.5
##   long NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1 20.5         0.033725924                 30751          0.23355016
## 2 20.5         0.015960753                 30751          0.16010043
## 3 20.5         0.010016616                 31111          0.12668130
## 4 20.5         0.001600598                 14563          0.05055509
## 5 20.5         0.018324556                  1958          0.17162744
## 6 20.5         0.004808260                 17177          0.08767895
##   numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## 1                 2         0.836709022                 31931
## 2                 2         0.894133247                   794
## 3                 2         0.661532267                   861
## 4                 2         0.265489647                 39602
## 5                 2         0.006257156                 39603
## 6                 2         0.036266868                 39426
##   NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today long_1800_today
## 1           1.3486467                 2           44.5           302.5
## 2           1.4077485                 2           44.5           302.5
## 3           1.1633920                 2           82.5           305.5
## 4           0.6854521                 2           53.5           172.5
## 5           0.1000496                 2           69.5            47.5
## 6           0.2423104                 2           61.5           188.5
##   latshift_1800_today latshift lat_today_2100_8.5 long_today_2100_8.5
## 1                  26       26               89.5               287.5
## 2                  27       27               89.5               288.5
## 3                 -10      -10               89.5               287.5
## 4                  20       20               85.5               279.5
## 5                   5        5               82.5               273.5
## 6                  14       14               82.5               271.5
# Visualize
Plot_nonInt(final_dat9$lat, final_dat9$long, 
            final_dat9$NN.sigma_today_2100_4.5, world, "sig nov 2100 4.5")

stationInfo9 <- stationInfo
names(stationInfo9) <- paste0(names(stationInfo),"_today_2100_4.5")
names(stationInfo9)[1] <- "NN.station_today_2100_4.5"
head(stationInfo9)
##   NN.station_today_2100_4.5 lat_today_2100_4.5 long_today_2100_4.5
## 1                         1              -69.5                20.5
## 2                         2              -68.5                20.5
## 3                         3              -67.5                20.5
## 4                         4              -66.5                20.5
## 5                         5              -65.5                20.5
## 6                         6              -64.5                20.5
final_dat10 <- left_join(final_dat9, stationInfo9)
## Joining, by = "NN.station_today_2100_4.5"
dim(final_dat10)
## [1] 39662    31
head(final_dat10)
##   No NN.sigma_today_2100_4.5 NN.station_today_2100_4.5 NN.Mdist_today_2100_4.5
## 1 36             0.260091777                     31342              0.67775158
## 2 37             0.150575943                     31338              0.50493740
## 3 38             0.005534165                     31456              0.09407847
## 4 39             0.139495959                     30347              0.48495881
## 5 40             0.011775140                      3437              0.13740012
## 6 41             0.098007656                      4227              0.40321232
##   numPCs_today_2100_4.5 NN.sigma_2100_8.5_today NN.station_2100_8.5_today
## 1                     2                8.292361                     40111
## 2                     2                8.292361                     40111
## 3                     2                8.292361                     40111
## 4                     2                8.292361                     40111
## 5                     2                8.292361                     40111
## 6                     2                8.292361                     40111
##   NN.Mdist_2100_8.5_today numPCs_2100_8.5_today NN.sigma_today_2100_8.5
## 1               10.055106                     2            0.3501681236
## 2               10.340310                     2            0.1043489359
## 3               10.400029                     2            0.0015820951
## 4               10.280321                     2            0.0006816715
## 5               10.091702                     2            0.1188184652
## 6                9.955604                     2            0.5463419465
##   NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5 numPCs_today_2100_8.5  lat
## 1                     29638              0.79989070                     2 70.5
## 2                     29695              0.41656957                     2 71.5
## 3                     29638              0.05026184                     2 72.5
## 4                     28978              0.03298615                     2 73.5
## 5                     28440              0.44577407                     2 74.5
## 6                     28254              1.03579199                     2 75.5
##   long NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 1 20.5         0.033725924                 30751          0.23355016
## 2 20.5         0.015960753                 30751          0.16010043
## 3 20.5         0.010016616                 31111          0.12668130
## 4 20.5         0.001600598                 14563          0.05055509
## 5 20.5         0.018324556                  1958          0.17162744
## 6 20.5         0.004808260                 17177          0.08767895
##   numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## 1                 2         0.836709022                 31931
## 2                 2         0.894133247                   794
## 3                 2         0.661532267                   861
## 4                 2         0.265489647                 39602
## 5                 2         0.006257156                 39603
## 6                 2         0.036266868                 39426
##   NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today long_1800_today
## 1           1.3486467                 2           44.5           302.5
## 2           1.4077485                 2           44.5           302.5
## 3           1.1633920                 2           82.5           305.5
## 4           0.6854521                 2           53.5           172.5
## 5           0.1000496                 2           69.5            47.5
## 6           0.2423104                 2           61.5           188.5
##   latshift_1800_today latshift lat_today_2100_8.5 long_today_2100_8.5
## 1                  26       26               89.5               287.5
## 2                  27       27               89.5               288.5
## 3                 -10      -10               89.5               287.5
## 4                  20       20               85.5               279.5
## 5                   5        5               82.5               273.5
## 6                  14       14               82.5               271.5
##   lat_today_2100_4.5 long_today_2100_4.5
## 1               86.5               307.5
## 2               82.5               307.5
## 3               85.5               308.5
## 4               69.5               298.5
## 5               70.5                61.5
## 6               72.5                68.5

today against 2100 RCP 4.5 (disappearing climate)

#--------------------------------  
### today against 2100 RCP 4.5 (disappearing climate) ####
#--------------------------------
NN.sigma_2100_4.5_today_N <- loop_sigma_D(B2[N_hem,], A1[N_hem,], C, "_2100_4.5_today")
## [1]   100   317 17413
## [1]   200   610 17413
## [1]   300   934 17413
## [1]   400  1247 17413
## [1]   500  1683 17413
## [1]   600  2053 17413
## [1]   700  2356 17413
## [1]   800  2660 17413
## [1]   900  2971 17413
## [1]  1000  3213 17413
## [1]  1100  3452 17413
## [1]  1200  3759 17413
## [1]  1300  3996 17413
## [1]  1400  4236 17413
## [1]  1500  4540 17413
## [1]  1600  4847 17413
## [1]  1700  5220 17413
## [1]  1800  5530 17413
## [1]  1900  5830 17413
## [1]  2000  6129 17413
## [1]  2100  6427 17413
## [1]  2200  6726 17413
## [1]  2300  7099 17413
## [1]  2400  7694 17413
## [1]  2500  8267 17413
## [1]  2600  8629 17413
## [1]  2700  8887 17413
## [1]  2800  9069 17413
## [1]  2900  9286 17413
## [1]  3000  9454 17413
## [1]  3100  9622 17413
## [1]  3200  9790 17413
## [1]  3300  9954 17413
## [1]  3400 10119 17413
## [1]  3500 10287 17413
## [1]  3600 10458 17413
## [1]  3700 10592 17413
## [1]  3800 10764 17413
## [1]  3900 10945 17413
## [1]  4000 11095 17413
## [1]  4100 11309 17413
## [1]  4200 11469 17413
## [1]  4300 11702 17413
## [1]  4400 11871 17413
## [1]  4500 12109 17413
## [1]  4600 12278 17413
## [1]  4700 12450 17413
## [1]  4800 12692 17413
## [1]  4900 12863 17413
## [1]  5000 13035 17413
## [1]  5100 13280 17413
## [1]  5200 13453 17413
## [1]  5300 13627 17413
## [1]  5400 13878 17413
## [1]  5500 14056 17413
## [1]  5600 14235 17413
## [1]  5700 14414 17413
## [1]  5800 14670 17413
## [1]  5900 14848 17413
## [1]  6000 15026 17413
## [1]  6100 15204 17413
## [1]  6200 15384 17413
## [1]  6300 15641 17413
## [1]  6400 15819 17413
## [1]  6500 15998 17413
## [1]  6600 16176 17413
## [1]  6700 16356 17413
## [1]  6800 16534 17413
## [1]  6900 16791 17413
## [1]  7000 16970 17413
## [1]  7100 17149 17413
## [1]  7200 17328 17413
## [1]  7300 17506 17413
## [1]  7400 17684 17413
## [1]  7500 17862 17413
## [1]  7600 18119 17413
## [1]  7700 18297 17413
## [1]  7800 18476 17413
## [1]  7900 18655 17413
## [1]  8000 18912 17413
## [1]  8100 19090 17413
## [1]  8200 19269 17413
## [1]  8300 19526 17413
## [1]  8400 19704 17413
## [1]  8500 19881 17413
## [1]  8600 20058 17413
## [1]  8700 20312 17413
## [1]  8800 20488 17413
## [1]  8900 20664 17413
## [1]  9000 20916 17413
## [1]  9100 21091 17413
## [1]  9200 21266 17413
## [1]  9300 21441 17413
## [1]  9400 21691 17413
## [1]  9500 21866 17413
## [1]  9600 22041 17413
## [1]  9700 22218 17413
## [1]  9800 22468 17413
## [1]  9900 22643 17413
## [1] 10000 22817 17413
## [1] 10100 22992 17413
## [1] 10200 23240 17413
## [1] 10300 23414 17413
## [1] 10400 23662 17413
## [1] 10500 23836 17413
## [1] 10600 24009 17413
## [1] 10700 24256 17413
## [1] 10800 24430 17413
## [1] 10900 24679 17413
## [1] 11000 24927 17413
## [1] 11100 25250 17413
## [1] 11200 25500 17413
## [1] 11300 25823 17413
## [1] 11400 26075 17413
## [1] 11500 26397 17413
## [1] 11600 26790 17413
## [1] 11700 27106 17413
## [1] 11800 27494 17413
## [1] 11900 27959 17413
## [1] 12000 28351 17413
## [1] 12100 28892 17413
## [1] 12200 29303 17413
## [1] 12300 29551 17413
## [1] 12400 29721 17413
## [1] 12500 29857 17413
## [1] 12600 30024 17413
## [1] 12700 30174 17413
## [1] 12800 30333 17413
## [1] 12900 30505 17413
## [1] 13000 30643 17413
## [1] 13100 30827 17413
## [1] 13200 30971 17413
## [1] 13300 31118 17413
## [1] 13400 31311 17413
## [1] 13500 31458 17413
## [1] 13600 31658 17413
## [1] 13700 31811 17413
## [1] 13800 32021 17413
## [1] 13900 32177 17413
## [1] 14000 32392 17413
## [1] 14100 32552 17413
## [1] 14200 32775 17413
## [1] 14300 32941 17413
## [1] 14400 33112 17413
## [1] 14500 33355 17413
## [1] 14600 33534 17413
## [1] 14700 33790 17413
## [1] 14800 33968 17413
## [1] 14900 34146 17413
## [1] 15000 34400 17413
## [1] 15100 34577 17413
## [1] 15200 34753 17413
## [1] 15300 35004 17413
## [1] 15400 35179 17413
## [1] 15500 35353 17413
## [1] 15600 35602 17413
## [1] 15700 35777 17413
## [1] 15800 35952 17413
## [1] 15900 36129 17413
## [1] 16000 36377 17413
## [1] 16100 36553 17413
## [1] 16200 36728 17413
## [1] 16300 36975 17413
## [1] 16400 37148 17413
## [1] 16500 37391 17413
## [1] 16600 37635 17413
## [1] 16700 37877 17413
## [1] 16800 38188 17413
## [1] 16900 38498 17413
## [1] 17000 38808 17413
## [1] 17100 39119 17413
## [1] 17200 39428 17413
## [1] 17300 39759 17413
## [1] 17400 40107 17413
NN.sigma_2100_4.5_today_S <- loop_sigma_D(B2[S_hem,], A1[S_hem,], C, "_2100_4.5_today")
## [1]   100   140 22249
## [1]   200   298 22249
## [1]   300   459 22249
## [1]   400   622 22249
## [1]   500   765 22249
## [1]   600   909 22249
## [1]   700  1052 22249
## [1]   800  1195 22249
## [1]   900  1345 22249
## [1]  1000  1468 22249
## [1]  1100  1619 22249
## [1]  1200  1763 22249
## [1]  1300  1923 22249
## [1]  1400  2093 22249
## [1]  1500  2228 22249
## [1]  1600  2401 22249
## [1]  1700  2539 22249
## [1]  1800  2718 22249
## [1]  1900  2858 22249
## [1]  2000  3040 22249
## [1]  2100  3182 22249
## [1]  2200  3372 22249
## [1]  2300  3518 22249
## [1]  2400  3710 22249
## [1]  2500  3856 22249
## [1]  2600  4047 22249
## [1]  2700  4190 22249
## [1]  2800  4368 22249
## [1]  2900  4506 22249
## [1]  3000  4678 22249
## [1]  3100  4810 22249
## [1]  3200  4967 22249
## [1]  3300  5092 22249
## [1]  3400  5240 22249
## [1]  3500  5367 22249
## [1]  3600  5531 22249
## [1]  3700  5666 22249
## [1]  3800  5834 22249
## [1]  3900  5969 22249
## [1]  4000  6141 22249
## [1]  4100  6278 22249
## [1]  4200  6451 22249
## [1]  4300  6586 22249
## [1]  4400  6753 22249
## [1]  4500  6882 22249
## [1]  4600  7031 22249
## [1]  4700  7149 22249
## [1]  4800  7275 22249
## [1]  4900  7388 22249
## [1]  5000  7514 22249
## [1]  5100  7640 22249
## [1]  5200  7753 22249
## [1]  5300  7882 22249
## [1]  5400  8010 22249
## [1]  5500  8139 22249
## [1]  5600  8256 22249
## [1]  5700  8404 22249
## [1]  5800  8538 22249
## [1]  5900  8680 22249
## [1]  6000  8864 22249
## [1]  6100  9073 22249
## [1]  6200  9259 22249
## [1]  6300  9501 22249
## [1]  6400  9758 22249
## [1]  6500 10019 22249
## [1]  6600 10336 22249
## [1]  6700 10614 22249
## [1]  6800 10899 22249
## [1]  6900 11128 22249
## [1]  7000 11363 22249
## [1]  7100 11532 22249
## [1]  7200 11769 22249
## [1]  7300 11938 22249
## [1]  7400 12178 22249
## [1]  7500 12350 22249
## [1]  7600 12523 22249
## [1]  7700 12772 22249
## [1]  7800 12948 22249
## [1]  7900 13203 22249
## [1]  8000 13382 22249
## [1]  8100 13561 22249
## [1]  8200 13824 22249
## [1]  8300 14012 22249
## [1]  8400 14279 22249
## [1]  8500 14462 22249
## [1]  8600 14649 22249
## [1]  8700 14915 22249
## [1]  8800 15102 22249
## [1]  8900 15289 22249
## [1]  9000 15556 22249
## [1]  9100 15740 22249
## [1]  9200 15926 22249
## [1]  9300 16112 22249
## [1]  9400 16386 22249
## [1]  9500 16573 22249
## [1]  9600 16760 22249
## [1]  9700 16947 22249
## [1]  9800 17224 22249
## [1]  9900 17414 22249
## [1] 10000 17604 22249
## [1] 10100 17881 22249
## [1] 10200 18065 22249
## [1] 10300 18247 22249
## [1] 10400 18428 22249
## [1] 10500 18688 22249
## [1] 10600 18867 22249
## [1] 10700 19044 22249
## [1] 10800 19298 22249
## [1] 10900 19475 22249
## [1] 11000 19652 22249
## [1] 11100 19828 22249
## [1] 11200 20081 22249
## [1] 11300 20258 22249
## [1] 11400 20436 22249
## [1] 11500 20694 22249
## [1] 11600 20874 22249
## [1] 11700 21054 22249
## [1] 11800 21314 22249
## [1] 11900 21494 22249
## [1] 12000 21674 22249
## [1] 12100 21934 22249
## [1] 12200 22115 22249
## [1] 12300 22295 22249
## [1] 12400 22553 22249
## [1] 12500 22731 22249
## [1] 12600 22909 22249
## [1] 12700 23161 22249
## [1] 12800 23335 22249
## [1] 12900 23508 22249
## [1] 13000 23753 22249
## [1] 13100 23923 22249
## [1] 13200 24092 22249
## [1] 13300 24326 22249
## [1] 13400 24485 22249
## [1] 13500 24690 22249
## [1] 13600 24839 22249
## [1] 13700 24987 22249
## [1] 13800 25178 22249
## [1] 13900 25322 22249
## [1] 14000 25465 22249
## [1] 14100 25646 22249
## [1] 14200 25785 22249
## [1] 14300 25923 22249
## [1] 14400 26095 22249
## [1] 14500 26228 22249
## [1] 14600 26360 22249
## [1] 14700 26522 22249
## [1] 14800 26652 22249
## [1] 14900 26811 22249
## [1] 15000 26938 22249
## [1] 15100 27065 22249
## [1] 15200 27218 22249
## [1] 15300 27344 22249
## [1] 15400 27470 22249
## [1] 15500 27619 22249
## [1] 15600 27743 22249
## [1] 15700 27890 22249
## [1] 15800 28011 22249
## [1] 15900 28132 22249
## [1] 16000 28274 22249
## [1] 16100 28394 22249
## [1] 16200 28514 22249
## [1] 16300 28648 22249
## [1] 16400 28764 22249
## [1] 16500 28898 22249
## [1] 16600 29022 22249
## [1] 16700 29172 22249
## [1] 16800 29327 22249
## [1] 16900 29498 22249
## [1] 17000 29708 22249
## [1] 17100 30060 22249
## [1] 17200 30380 22249
## [1] 17300 30680 22249
## [1] 17400 30923 22249
## [1] 17500 31234 22249
## [1] 17600 31474 22249
## [1] 17700 31716 22249
## [1] 17800 31954 22249
## [1] 17900 32192 22249
## [1] 18000 32359 22249
## [1] 18100 32598 22249
## [1] 18200 32844 22249
## [1] 18300 33019 22249
## [1] 18400 33266 22249
## [1] 18500 33439 22249
## [1] 18600 33616 22249
## [1] 18700 33863 22249
## [1] 18800 34038 22249
## [1] 18900 34213 22249
## [1] 19000 34463 22249
## [1] 19100 34638 22249
## [1] 19200 34814 22249
## [1] 19300 35066 22249
## [1] 19400 35242 22249
## [1] 19500 35419 22249
## [1] 19600 35678 22249
## [1] 19700 35859 22249
## [1] 19800 36042 22249
## [1] 19900 36312 22249
## [1] 20000 36490 22249
## [1] 20100 36664 22249
## [1] 20200 36908 22249
## [1] 20300 37078 22249
## [1] 20400 37310 22249
## [1] 20500 37464 22249
## [1] 20600 37612 22249
## [1] 20700 37801 22249
## [1] 20800 37943 22249
## [1] 20900 38122 22249
## [1] 21000 38258 22249
## [1] 21100 38428 22249
## [1] 21200 38563 22249
## [1] 21300 38697 22249
## [1] 21400 38864 22249
## [1] 21500 38997 22249
## [1] 21600 39161 22249
## [1] 21700 39292 22249
## [1] 21800 39446 22249
## [1] 21900 39572 22249
## [1] 22000 39721 22249
## [1] 22100 39866 22249
## [1] 22200 40030 22249
NN.sigma_2100_4.5_today <- rbind(NN.sigma_2100_4.5_today_N, NN.sigma_2100_4.5_today_S)


which(is.infinite(NN.sigma_2100_4.5_today$NN.sigma_2100_4.5_today))
## integer(0)
which(is.na(NN.sigma_2100_4.5_today$NN.sigma_2100_4.5_today))
## integer(0)
final_dat11 <- full_join(NN.sigma_2100_4.5_today, final_dat10, by="No")
head(final_dat11)
##   No NN.sigma_2100_4.5_today NN.station_2100_4.5_today NN.Mdist_2100_4.5_today
## 1 36                3.551980                     40055                3.967162
## 2 37                3.505984                     40055                3.923107
## 3 38                3.315848                     40055                3.741139
## 4 39                2.951652                     40055                3.393209
## 5 40                2.553620                     40111                3.013691
## 6 41                2.223800                     40111                2.699426
##   numPCs_2100_4.5_today NN.sigma_today_2100_4.5 NN.station_today_2100_4.5
## 1                     2             0.260091777                     31342
## 2                     2             0.150575943                     31338
## 3                     2             0.005534165                     31456
## 4                     2             0.139495959                     30347
## 5                     2             0.011775140                      3437
## 6                     2             0.098007656                      4227
##   NN.Mdist_today_2100_4.5 numPCs_today_2100_4.5 NN.sigma_2100_8.5_today
## 1              0.67775158                     2                8.292361
## 2              0.50493740                     2                8.292361
## 3              0.09407847                     2                8.292361
## 4              0.48495881                     2                8.292361
## 5              0.13740012                     2                8.292361
## 6              0.40321232                     2                8.292361
##   NN.station_2100_8.5_today NN.Mdist_2100_8.5_today numPCs_2100_8.5_today
## 1                     40111               10.055106                     2
## 2                     40111               10.340310                     2
## 3                     40111               10.400029                     2
## 4                     40111               10.280321                     2
## 5                     40111               10.091702                     2
## 6                     40111                9.955604                     2
##   NN.sigma_today_2100_8.5 NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5
## 1            0.3501681236                     29638              0.79989070
## 2            0.1043489359                     29695              0.41656957
## 3            0.0015820951                     29638              0.05026184
## 4            0.0006816715                     28978              0.03298615
## 5            0.1188184652                     28440              0.44577407
## 6            0.5463419465                     28254              1.03579199
##   numPCs_today_2100_8.5  lat long NN.sigma_1800_today NN.station_1800_today
## 1                     2 70.5 20.5         0.033725924                 30751
## 2                     2 71.5 20.5         0.015960753                 30751
## 3                     2 72.5 20.5         0.010016616                 31111
## 4                     2 73.5 20.5         0.001600598                 14563
## 5                     2 74.5 20.5         0.018324556                  1958
## 6                     2 75.5 20.5         0.004808260                 17177
##   NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 1          0.23355016                 2         0.836709022
## 2          0.16010043                 2         0.894133247
## 3          0.12668130                 2         0.661532267
## 4          0.05055509                 2         0.265489647
## 5          0.17162744                 2         0.006257156
## 6          0.08767895                 2         0.036266868
##   NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## 1                 31931           1.3486467                 2           44.5
## 2                   794           1.4077485                 2           44.5
## 3                   861           1.1633920                 2           82.5
## 4                 39602           0.6854521                 2           53.5
## 5                 39603           0.1000496                 2           69.5
## 6                 39426           0.2423104                 2           61.5
##   long_1800_today latshift_1800_today latshift lat_today_2100_8.5
## 1           302.5                  26       26               89.5
## 2           302.5                  27       27               89.5
## 3           305.5                 -10      -10               89.5
## 4           172.5                  20       20               85.5
## 5            47.5                   5        5               82.5
## 6           188.5                  14       14               82.5
##   long_today_2100_8.5 lat_today_2100_4.5 long_today_2100_4.5
## 1               287.5               86.5               307.5
## 2               288.5               82.5               307.5
## 3               287.5               85.5               308.5
## 4               279.5               69.5               298.5
## 5               273.5               70.5                61.5
## 6               271.5               72.5                68.5
# Visualize
Plot_nonInt(final_dat11$lat, final_dat11$long, 
            final_dat11$NN.sigma_2100_4.5_today, world, "sigma dis. 2000 in 2100 4.5")

# stationInfo11 <- stationInfo
# names(stationInfo11) <- paste0(names(stationInfo),"_2100_4.5_today")
# names(stationInfo11)[1] <- "NN.station_2100_4.5_today"
# head(stationInfo11)
# final_dat12 <- left_join(final_dat11, stationInfo11)
final_dat12 <- final_dat11
dim(final_dat12)
## [1] 39662    35
head(final_dat12)
##   No NN.sigma_2100_4.5_today NN.station_2100_4.5_today NN.Mdist_2100_4.5_today
## 1 36                3.551980                     40055                3.967162
## 2 37                3.505984                     40055                3.923107
## 3 38                3.315848                     40055                3.741139
## 4 39                2.951652                     40055                3.393209
## 5 40                2.553620                     40111                3.013691
## 6 41                2.223800                     40111                2.699426
##   numPCs_2100_4.5_today NN.sigma_today_2100_4.5 NN.station_today_2100_4.5
## 1                     2             0.260091777                     31342
## 2                     2             0.150575943                     31338
## 3                     2             0.005534165                     31456
## 4                     2             0.139495959                     30347
## 5                     2             0.011775140                      3437
## 6                     2             0.098007656                      4227
##   NN.Mdist_today_2100_4.5 numPCs_today_2100_4.5 NN.sigma_2100_8.5_today
## 1              0.67775158                     2                8.292361
## 2              0.50493740                     2                8.292361
## 3              0.09407847                     2                8.292361
## 4              0.48495881                     2                8.292361
## 5              0.13740012                     2                8.292361
## 6              0.40321232                     2                8.292361
##   NN.station_2100_8.5_today NN.Mdist_2100_8.5_today numPCs_2100_8.5_today
## 1                     40111               10.055106                     2
## 2                     40111               10.340310                     2
## 3                     40111               10.400029                     2
## 4                     40111               10.280321                     2
## 5                     40111               10.091702                     2
## 6                     40111                9.955604                     2
##   NN.sigma_today_2100_8.5 NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5
## 1            0.3501681236                     29638              0.79989070
## 2            0.1043489359                     29695              0.41656957
## 3            0.0015820951                     29638              0.05026184
## 4            0.0006816715                     28978              0.03298615
## 5            0.1188184652                     28440              0.44577407
## 6            0.5463419465                     28254              1.03579199
##   numPCs_today_2100_8.5  lat long NN.sigma_1800_today NN.station_1800_today
## 1                     2 70.5 20.5         0.033725924                 30751
## 2                     2 71.5 20.5         0.015960753                 30751
## 3                     2 72.5 20.5         0.010016616                 31111
## 4                     2 73.5 20.5         0.001600598                 14563
## 5                     2 74.5 20.5         0.018324556                  1958
## 6                     2 75.5 20.5         0.004808260                 17177
##   NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 1          0.23355016                 2         0.836709022
## 2          0.16010043                 2         0.894133247
## 3          0.12668130                 2         0.661532267
## 4          0.05055509                 2         0.265489647
## 5          0.17162744                 2         0.006257156
## 6          0.08767895                 2         0.036266868
##   NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## 1                 31931           1.3486467                 2           44.5
## 2                   794           1.4077485                 2           44.5
## 3                   861           1.1633920                 2           82.5
## 4                 39602           0.6854521                 2           53.5
## 5                 39603           0.1000496                 2           69.5
## 6                 39426           0.2423104                 2           61.5
##   long_1800_today latshift_1800_today latshift lat_today_2100_8.5
## 1           302.5                  26       26               89.5
## 2           302.5                  27       27               89.5
## 3           305.5                 -10      -10               89.5
## 4           172.5                  20       20               85.5
## 5            47.5                   5        5               82.5
## 6           188.5                  14       14               82.5
##   long_today_2100_8.5 lat_today_2100_4.5 long_today_2100_4.5
## 1               287.5               86.5               307.5
## 2               288.5               82.5               307.5
## 3               287.5               85.5               308.5
## 4               279.5               69.5               298.5
## 5               273.5               70.5                61.5
## 6               271.5               72.5                68.5
sum(!complete.cases(final_dat12))
## [1] 22249
### Final dataframe ####
dat_Nov <- final_dat12
### Write to file ####
write.csv(final_dat12, paste0(results_dir, "SigmaD_Hem.csv"), row.names = FALSE)

head(dat_Nov)
##   No NN.sigma_2100_4.5_today NN.station_2100_4.5_today NN.Mdist_2100_4.5_today
## 1 36                3.551980                     40055                3.967162
## 2 37                3.505984                     40055                3.923107
## 3 38                3.315848                     40055                3.741139
## 4 39                2.951652                     40055                3.393209
## 5 40                2.553620                     40111                3.013691
## 6 41                2.223800                     40111                2.699426
##   numPCs_2100_4.5_today NN.sigma_today_2100_4.5 NN.station_today_2100_4.5
## 1                     2             0.260091777                     31342
## 2                     2             0.150575943                     31338
## 3                     2             0.005534165                     31456
## 4                     2             0.139495959                     30347
## 5                     2             0.011775140                      3437
## 6                     2             0.098007656                      4227
##   NN.Mdist_today_2100_4.5 numPCs_today_2100_4.5 NN.sigma_2100_8.5_today
## 1              0.67775158                     2                8.292361
## 2              0.50493740                     2                8.292361
## 3              0.09407847                     2                8.292361
## 4              0.48495881                     2                8.292361
## 5              0.13740012                     2                8.292361
## 6              0.40321232                     2                8.292361
##   NN.station_2100_8.5_today NN.Mdist_2100_8.5_today numPCs_2100_8.5_today
## 1                     40111               10.055106                     2
## 2                     40111               10.340310                     2
## 3                     40111               10.400029                     2
## 4                     40111               10.280321                     2
## 5                     40111               10.091702                     2
## 6                     40111                9.955604                     2
##   NN.sigma_today_2100_8.5 NN.station_today_2100_8.5 NN.Mdist_today_2100_8.5
## 1            0.3501681236                     29638              0.79989070
## 2            0.1043489359                     29695              0.41656957
## 3            0.0015820951                     29638              0.05026184
## 4            0.0006816715                     28978              0.03298615
## 5            0.1188184652                     28440              0.44577407
## 6            0.5463419465                     28254              1.03579199
##   numPCs_today_2100_8.5  lat long NN.sigma_1800_today NN.station_1800_today
## 1                     2 70.5 20.5         0.033725924                 30751
## 2                     2 71.5 20.5         0.015960753                 30751
## 3                     2 72.5 20.5         0.010016616                 31111
## 4                     2 73.5 20.5         0.001600598                 14563
## 5                     2 74.5 20.5         0.018324556                  1958
## 6                     2 75.5 20.5         0.004808260                 17177
##   NN.Mdist_1800_today numPCs_1800_today NN.sigma_today_1800
## 1          0.23355016                 2         0.836709022
## 2          0.16010043                 2         0.894133247
## 3          0.12668130                 2         0.661532267
## 4          0.05055509                 2         0.265489647
## 5          0.17162744                 2         0.006257156
## 6          0.08767895                 2         0.036266868
##   NN.station_today_1800 NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today
## 1                 31931           1.3486467                 2           44.5
## 2                   794           1.4077485                 2           44.5
## 3                   861           1.1633920                 2           82.5
## 4                 39602           0.6854521                 2           53.5
## 5                 39603           0.1000496                 2           69.5
## 6                 39426           0.2423104                 2           61.5
##   long_1800_today latshift_1800_today latshift lat_today_2100_8.5
## 1           302.5                  26       26               89.5
## 2           302.5                  27       27               89.5
## 3           305.5                 -10      -10               89.5
## 4           172.5                  20       20               85.5
## 5            47.5                   5        5               82.5
## 6           188.5                  14       14               82.5
##   long_today_2100_8.5 lat_today_2100_4.5 long_today_2100_4.5
## 1               287.5               86.5               307.5
## 2               288.5               82.5               307.5
## 3               287.5               85.5               308.5
## 4               279.5               69.5               298.5
## 5               273.5               70.5                61.5
## 6               271.5               72.5                68.5
for_liqing <- data.frame(No=dat_Nov$No, lat=dat_Nov$lat,long=dat_Nov$long,
                         A=dat_Nov$NN.sigma_today_1800, 
                         B=dat_Nov$NN.sigma_1800_today, 
                         C=dat_Nov$NN.sigma_2100_4.5_today, 
                         D=dat_Nov$NN.sigma_today_2100_4.5, 
                         E=dat_Nov$NN.sigma_2100_8.5_today, 
                         F=dat_Nov$NN.sigma_today_2100_8.5)
write.csv(for_liqing,  paste0(results_dir,"SigmaD_Hem_plot1.csv"), row.names=FALSE)

for_liqing_Md <- data.frame(No=dat_Nov$No, lat=dat_Nov$lat,long=dat_Nov$long,A=dat_Nov$NN.Mdist_today_1800, B=dat_Nov$NN.Mdist_1800_today, C=dat_Nov$NN.Mdist_2100_4.5_today, D=dat_Nov$NN.Mdist_today_2100_4.5, E=dat_Nov$NN.Mdist_2100_8.5_today, F=dat_Nov$NN.Mdist_today_2100_8.5)
write.csv(for_liqing_Md,  paste0(results_dir,"SigmaD_Hem_plot2.csv"), row.names=FALSE)
#plot(final_dat12$lat_today_2100_8.5, final_dat12$lat) #If A is today and B is future, 
#  where station No's climate in the future will come from today (or the closest similar climate today).
# Latitude of the station today (x-axis) where the latitude of the station in the future (y-axis) will come from


#plot(final_dat12$lat, final_dat12$lat_2100_8.5_today) # If A is future and B is today, this 
# represents where station "No" will be found in the future (or the closest similar climate).
# Latitude of the station today (x-axis) and where it will be in the future (y-axis)

Calculate percentages

### Calculate percentages ####  
N <- nrow(dat_Nov)  

calc_perc <- function(x, descrip){
  print(paste("Percent of cells with moderate dissimilarity for", descrip,":", 
              round(sum(x>2 & x<4)/N,3)*100, sep=" "))
  print(paste("Percent of cells with extreme dissimilarity for", descrip,":", 
              round(sum(x>4)/N, 3)*100, sep=" "))
}

calc_perc(dat_Nov$NN.sigma_today_1800, "1800 disappearing in 2000")
## [1] "Percent of cells with moderate dissimilarity for 1800 disappearing in 2000 : 12.4"
## [1] "Percent of cells with extreme dissimilarity for 1800 disappearing in 2000 : 0"
calc_perc(dat_Nov$NN.sigma_1800_today, "Novel climates in 2000 since 1800")
## [1] "Percent of cells with moderate dissimilarity for Novel climates in 2000 since 1800 : 3.7"
## [1] "Percent of cells with extreme dissimilarity for Novel climates in 2000 since 1800 : 0"
calc_perc(dat_Nov$NN.sigma_2100_4.5_today, "2000 disappearing in 2100 4.5")
## [1] "Percent of cells with moderate dissimilarity for 2000 disappearing in 2100 4.5 : 34.2"
## [1] "Percent of cells with extreme dissimilarity for 2000 disappearing in 2100 4.5 : 35.6"
calc_perc(dat_Nov$NN.sigma_today_2100_4.5, "Novel climates in 2100 4.5 since 2000")
## [1] "Percent of cells with moderate dissimilarity for Novel climates in 2100 4.5 since 2000 : 42.1"
## [1] "Percent of cells with extreme dissimilarity for Novel climates in 2100 4.5 since 2000 : 10.3"
calc_perc(dat_Nov$NN.sigma_2100_8.5_today, "2000 disappearing in 2100 8.5")
## [1] "Percent of cells with moderate dissimilarity for 2000 disappearing in 2100 8.5 : 3.5"
## [1] "Percent of cells with extreme dissimilarity for 2000 disappearing in 2100 8.5 : 95"
calc_perc(dat_Nov$NN.sigma_today_2100_8.5, "Novel climates in 2100 8.5 since 2000")
## [1] "Percent of cells with moderate dissimilarity for Novel climates in 2100 8.5 since 2000 : 6.7"
## [1] "Percent of cells with extreme dissimilarity for Novel climates in 2100 8.5 since 2000 : 81.9"
# Visualize relationship between sigma_D and M_D
plot(dat_Nov$NN.sigma_today_1800, dat_Nov$NN.Mdist_today_1800)

plot(dat_Nov$NN.sigma_today_2100_8.5, dat_Nov$NN.Mdist_today_2100_8.5, col=dat_Nov$numPCs_today_2100_8.5)
legend(0, 15, 1:4, 1:4)
abline(h=12)

hist(dat_Nov$numPCs_today_2100_8.5)